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ABSTRACT 


The  turbulent  free  v  .tvection  of  air  above  a  2-foot  diameter, 
heated  horizontal  plate  has  been  studied  experimentally  and  numerically. 
The  mean  temperature  fields  and  the  indraft  profiles  for  two  mean  plate 
temperatures  were  measured  using  a  thermocouple  and  a  constant  tempera¬ 
ture  hot-wire  anemometer.  Also,  the  turbulence  and  mean  velocity  were 
measured  for  the  higher  plate  temperature  using  the  hot-wire  method. 

The  flow  field  was  visualized  by  shadow  photograph  technique.  From 
visualization  and  measurements,  it  was  found  that  the  region  of  signi¬ 
ficant  deviation  from  ambient  t  '.mperature  and  velocity  was  restricted 
to  a  region  near  the  plate  centerline  (the  primary  flow  region).  The 
indraft  velocity  was  found  to  be  relatively  large  near  the  ground  level 
(within  approximately  1"  of  the  ground). 

The  major  temperature  drop  took  place  in  the  region  very  near  the 
plate.  Within  0.02"  of  the  plate  the  temperature  distribution  in  the 

air  could  be  calculated  based  on  conduction  only.  This  region  was  _ 

therefore,  called  the  "conduction  layer."  At  a  given  mean  plate  tempera¬ 
ture,  the  temperature  gradient  was  found  to  increase  with  the  radius. 

Data  obtained  from  heat-transfer  measurements  were  consistent  with  the 
one-third  power  correlation  reported  in. the  literature. 

The  turbulence  in  the  flow  field  was  found  to  consist  of  low 
frequency  and  high  amplitude  fluctuations  (on  the  order  of  10  Hz  and 
1  ft/sec).  Because  of  the  limitation  of  the  hot-wire  technique  for 
large  turbulence  measurements,  flow  velocities  could  not  be  deduced 
directly  from  hot-wire  data.  To  remove  this  difficulty,  a  numerical 
data  simulation  scheme  has  been  developed  in  which  the  parameters 
describing  the  turbulent  flow  (r.m.s.  fluctuations  and  correlation 
coefficients)  were  used  as  input.  By  inferring  from  the  simulated  data 
of  know  parameters,  experimental  hot-wire  data  reduction  was  then 
possible.  Data  reduction  model  was  validated  by  numerical  experiments. 

The  eddy  diffusivity  in  the  region  away  from  the  conduction  layer 
was  estimated  based  on  temperature,  velocity  and  turbulence  data  using 
two  independent  methods.  The  agreement  was  good.  The  spatial  varia¬ 
tions  of  the  eddy  diffusivity  in  most  of  the  primary  flow  region  was 


found  to  be  gradual  with  rapid  drops  occurring  in  the  region  between 
the  primary  flow  and  the  cold  ambient. 

A  numerical  flow  calculation  was  made.  The  mathematical  formula¬ 
tion  was  baped  on  Boussinesq  approximations  using  a  constant  eddy 
diffusivity  model.  A  turbulent  Grashof  number  Gr-j  (the  governing 
parameter)  was  defined  through  the  definition  of  a  characteristic  plate 
temperature  rise  the  plate  mean  heat  flux  and  the  eddy  diffusivity. 

Gr^  and  were  obtained  based  on  the  best  fit  of  experimental  and 
numerical  centerline  temperatures. 

By  the  specification  of  at  the  plate  surface,  the  effect  of 
the  intense  variation  of  eddy  diffusivity  in  the  conduction  layer  region 
could  be  avoided  in  the  numerical  calculations.  Numerical  results 
based  on  a  constant  eddy  diffusivity  model  were  obtained  and  compared 
with  the  experimental  data.  Due  apparently  to  the  non-constancy  of 
the  eddy  diffusivity,  the  calculated  temperature  and  velocity  profiles 
exhibit  less  constriction  than  the  experimental  data.  Therefore  a  more 
general  turbulent  transport  model  will  be  required  to  provide  a  good 
theoretical  description  of  the  phenomena. 
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I.  INTRODUCTION 


Theoretical  investigations  of  free  convection  problems  have  been 
concerned  primarily  with  problems  of  simple  geometry  such  at  a  fluid 
layer  above  a  large,  horizontal,  heated  surface,  or  one  located  between 
two  large,  horizontal,  parallel  surfaces.  These  geometries  allow  one- 
dimensional  formulations.  Using  Boussiaesq  approximations  (1), 

Howard  (2)  obtained  a  solution  for  the  one-dimensional  turbulent  convex 
tion  above  a  heated  horizontal  surface.  He  has  demonstrated  agreement 
between  his  solution  end  the  experimental  data  of  Townsend  (3). 

At  some  distance  above  the  surface  where  the  updraft  is  a  dominant 
feature,  the  flow  resembles  a  free  jet.  This  flow  region  is  called  the 
convection  plume.  The  boundary  layer  approximations  are  then  applicable 
(4).  Similarity  solutions  have  been  obtained  for  problems  such  as 
convection  above  a  line  heat  source  or  a  linear  heat  source  of  finite 
width.  The  results  are  in  reasonable  agreement  with  experimental 
measurements  (e.g.,  reference  5,  6  and  7).  Because  of  the  conditions 
of  zero  average  shear  stress  and  heat  flux  along  the  plume  axis  and  at 
infinity,  the  shear  stress  and  heat  flux  terms  appearing  in  the  governing 
differential  equations  may  be  eliminated  by  integration  in  the  trans¬ 
verse  direction  from  the  plume  axis  to  infinity.  The  solution  to  the 
problem  may  be  obtained  by  use  of  either  assumed  or  measured  profile 
functions  of  mean  temperature  and  velocity  in  the  transverse  direction. 


Numbers  in  parentheses  indicate  references  listed  in 
the  Bibliography. 


Since  nc-  explicit  knowledge  of  the  turbulence  structure  is  necessary 
for  obtaining  the  plume  solutions,  the  study  of  turbulence  structure 
even  in  the  plume  region  has  been  only  in  the  infant  stage  (8) . 

For  multi-dimensional  flows  of  the  type  studied  in  this  investiga¬ 
tion,  the  highly  non-linear  nature  of  the  governing  equations  effectively 
precludes  analytical  solutions.  Recent  development  in  numerical  methods 
utilizing  large  digital  computers  has  made  possible  the  treatment  of 
some  of  the  non-linear  features  in  the  equations  of  motion,  for  such 
problems  as  convection  in  an  enclosed  cavity  (9) ,  and  wake  flow  behind 
an  obstacle  (10).  Reasonably  successful  solutions  can  usually  be 
obtained  for  problems  of  this  nature.  An  additional  uncertainty  in  the 
problem  formulation  for  the  work  described  here  derives  from  the  fact 
that  an  infinite  domain  must  be  approximated  by  a  finite  one  and  it  is 
not  clear  how  to  specify  boundary  conditions  at  the  artificial  boundaries. 
Beyond  this,  it  is  in  general  not  easy  to  demonstrate  that  the  results 

ui  iiutaci'iCax  caicuxai>iOuS  Of  COuiplCX  nOtl*'liu£a&  phcriviTictla  actually 

constitute  a  valid  solution  of  the  problem  as  formulated. 

In  view  of  these  multiple  uncertainties,  the  proper  initial 
approach  is  to  build  understanding  on  the  basis  of  a  simple  numerical 
calculation  model  which  is  easily  interpreted  in  the  light  of  experi¬ 
mental  evidence  and  which  can  be  reduced  to  analytically  tractable 
limiting  cases  for  validating  the  calculations.  This  approach  has, 
to  varying  degrees,  been  followed  for  several  free  convection  problems 
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(e.g.,  references  11,  12,  13).  Among  these,  only  one  (13)  describes 
numerical  flow  field  calculations  directly  relevant  to  the  present 
problem.  Unfortunately,  no  correlation  with  experimental  data  was 
reported. 

Although  there  exists  a  large  body  of  experimental  literature 
concerning  free  convection  in  general,  very  little  is  reported  which  is 
directly  applicable  to  the  present  study.  The  particularly  difficult 
experimental  feature  of  this  problem  is  the  accurate  measurement  of 
the  flow  velocity  vector.  For  any  reasonable  experimental  dimensions, 
velocities  near  the  heater  plate  are  small,  i.e.,  on  the  order  of  a 
fraction  of  a  foot  per  second.  Moreover,  turbulent  fluctuation 
velocities  in  some  regions  of  the  flow  field  are  of  the  same  order  of 
magnitude  as  the  mean.  In  the  experimental  portion  of  this  investigation, 
a  hot-wire  anemometer  was  used  to  obtain  flow  field  data.  Toward  the 
low  end  of  the  velocity  range  of  interest,  hot-wire  techniques  suffer 
from  orientation  effects  due  to  local  perturbation  of  the  flow  by  the 
hot-wire  itself.  Recently,  laser  doppler  techniques  (14,15)  have  been 
developed  which  may  obviate  most  of  the  difficulties  associated  with 
hot-wire  techniques.  However,  these  techniques  appeared  after  the  hot¬ 
wire  measurements  reported  in  this  thesis  were  well  underway,  and  were 
not  used  for  the  present  work. 

SCOPS  OF  THIS  WORK 

This  thesis  consists  of  two  major  parts:  experimental  measurements 
and  numerical  flow  and  temperature  field  calculations.  A  constant 
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temperature  hot-wire  anemometer  together  with  a  thermocouple  was  used 
for  obtaining  the  velocity  and  temperature  data.  The  numerical  calcula¬ 
tion  was  intended  for  exploring  the  qualitative  features  of  the  flow 
field.  Therefore,  calculations  were  made  based  on  the  Boussineaq 
approximations  (1)  incorporating  a  constant  eddy  diffusivity  model. 

With  the  mesh  spacing  (0.75  inch)  used  in  the  study,  the  flow  between 
the  ground  and  the  nodel  points  next  to  it  cannot  be  predicted.  Since 
the  indraft  at  the  ground  level  is  one  of  the  most  interesting  features 
in  the  flow  field,  the  numerical  results  must  be  supplemented  to  make 
possible  the  indraft  prediction  very  near  the  ground.  This  was 
accomplished  by  an  integral  method. 


II.  LITERATURE  SURVEY 


The  quantitative  investigation  of  a  fluid  layer  heated  from  below 
may  be  traced  back  to  the  late  1800's  when  Thomson  (16)  and  later 
Benard  (17)  observed  the  flow  patterns  following  the  onset  of  thermal 
instability.  Lord  Rayleigh  (18)  used  a  linearized  stability  theory  and 
established  the  criterion  (a  critical  Rayleigh  Number)  at  which  convective 
motion  commences.  Since  then  the  linear  instability  properties  for 
fluid  layers  of  infinite  horizontal  extent  have  been  studied  by  a 
number  of  investigators.  A  comprehensive  account  of  this  type  of 
problem  appears  in  Chandrasekhar's  book  (19). 

The  convection  which  develops  due  to  a  heated  surface  of  finite 
size  is  the  subject  of  concern  in  this  investigation.  Due  to  the  three- 
dimensionality  of  the  flow  configuration,  this  problem  is  quite  different 
from  that  due  to  an  infinite  surface.  Linearization  has  not  been 
possible,  and  no  analytical  approach  to  this  problem  has  been  reported 
in  literature. 

Examination  of  the  available  literature  reveals  that  knowledge  of 
the  turbulent  flow  field  above  a  heated  horizontal  plate  of  finite 
extent  exists  only  in  a  primitive  state.  This  may  be  attributed  to  the 
intrinsic  non-linearity  of  the  basic  governing  equations.  It  is 
complicated  further  by  the  lack  of  an  adequate  description  of  the 
turbulent  transport  mechanism.  Relevant  experimental,  mathematical 


and  numerical  studies  are  discussed  below. 


EXPERIMENTAL  STUDIES 


The  measurement  of  free  convection  above  a  heated  horizontal  surface 
was  apparently  first  reported  by  Ramdas  (20,21).  He  measured  temperature 
variations  above  heated  land  surface  (e.g.,  during  the  day)  and  concluded 
that  the  variation  of  temperature  was  most  rapid  nearest  to  the  surface, 
and  was  linear.  Information  of  this  type  is  important  to  the  under¬ 
standing  of  mass  movement  of  the  atmospheric  air  (vind  formation)  and 
many  early  investigations  were  made  by  meteorologists  (e.g.,  references 
22,  23,  24). 

Heat  transfer  data  of  an  engineering  nature  were  presented  by 
Schmidt  (25),  who  photographed  the  boundary  layer  formed  in  the 
neighborhood  of  a  heated  horizontal  plate.  Later,  Fishenden  and 
Saunders  (26)  measured  the  heat  transfer  from  horizontal,  rectangular 
plates  of  sizes  up  to  2-foot  square  and  at  temperature  up  to  1000°F  above 
the  ambient.  They  presented  data  for  both  the  cases  of  a  heated  surface 
facing  up  and  facing  down  but  gave  no  detailed  description  of  the 
experimental  configuration.  The  data  were  correlated  in  the  familiar 
fashion  of  Nusselt  Number  (Nu)  versus  Rayleigh  Number  (Ra) .  For  the 

turbulent  range,  the  data  gave  the  correlation: 

1/3 

Nu  =  constant  (Ra)  ' 

Since  Nu  contains  a  length  scale  and  the  heat  flux  to  the  first 
power  and  Ra  contains  a  length  scale  to  the  third  power,  the  above 
correlation  implies  that  the  heat  flux  is  independent  of  length  scale 
in  the  turbulent  regime.  A  comprehensive  summary  of  the  heat  transfer 
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data  from  horizontal  plate  was  compiled  by  Jakob  (27),  however,  little 
in  the  way  of  experimental  detail  and  measurement  technique  were 
discussed. 

Townsend  (3)  measured  the  temperature  fluctuations  above  the 
middle  of  the  uniformly  heated  bottom  of  an  open-topped  box  (30cc  x  40ca 
x  56cm  high).  All  the  measurements  were  made  within  8  cm  of  the  bottom 
surface.  The  mean  temperature  distribution  near  the  surface  was 
reasonably  linear.  From  the  implied  temperature  gradient  and  the  overall 
heat  transfer,  he  obtained  the  same  correlation  as  did  Fishenden  and 
Saunders  (26).  The  data  showed  relatively  large  scatter  when  the  surface 
temperature  was  high.  The  largest  temperature  fluctuations  occurred 
within  0.5  cm  of  the  surface  and  were  about  20%  of  the  temperature 
difference  between  the  plate  and  the  ambient.  The  rate  of  spatial 
temperature  decay  was  found  to  decrease  with  plate  temperature. 

Tritton  (28)  investigated  the  turbulent  free  convection  above  a 
heated  plate  inclined  at  a  small  angle  to  the  horizontal.  He  used  a 
resistance  wire  thermometer  to  measure  the  temperature  and  a  quartz 
fiber  anemometer  to  measure  the  velocity  parallel  to  the  plate  in  the 
boundary  layer  region.  He  found  that  the  temperature  field  was  not 
greatly  altered  due  to  place  inclination  (no  systematic  variations  were 
observed)  and  that  the  mean  temperature  field  was  largely  controlled  by 
the  turbulence.  He  succeeded  in  measuring  the  velocity  in  the  boundary 
layer  region  by  observing  the  deflection  of  a  cantilevered  quartz  fiber 
in  the  flow  field.  This  method,  however,  was  apparently  handicapped  by 
the  lack  of  a  continuous  recording  technique. 
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Project  FLAMBEAU  (29)  is  a  program  of  research  on  large  area  fires. 
During  recent  years,  a  large  amount  of  experimental  data  (e.g.,  fuel 
consumption  rate,  flow  and  temperature  field  information,  radiation 
effects,  etc.)  has  been  collected  from  large  size  experimental  fires  of 
horizontal  length  scale  up  to  1000  feet.  Little  systematic  analysis  cf 
these  data  has  been  reported. 

Recently,  Parker,  Corlett  and  Lee  (30)  attempted  to  test  the 
hypothesis  that  Grashof  number  is  unimportant  in  determining  such  flow 
field  characteristics  as  inflow  velocities.  Using  a  hot-wire  anemometer, 
they  measured  the  inflow  velocity  at  the  edge  of  a  square  array  (24"  x 
24")  of  horizontally  arranged  electric  heater  elements  which,  according 
to  the  hypothesis,  should  scale  the  September  1967  FLAMBEAU  test  fire. 
Insofar  as  the  severely  limited  prototype  data  permitted  comparison,  the 
hypothesis  was  found  correct.  Lore  data  of  the  same  work  in  improved 
form  may  be  found  in  reference  31. 

MATHEMATICAL  STUDIES 

Boussinesq  (1)  in  1903  introduced  nis  well-known  approximations  for 
calculating  thermal  convection  problems.  The  Boussinesq  approximations 
are:  (1)  Density  variation  due  to  pressure  (as  opposed  to  temperature) 
variation  is  legligible;  (2)  Density  variation  in  the  governing 
differential  equations  is  negligible  except  in  the  buoyancy  term  in  the 
momentum  equation. 

The  mathematical  formulation  of  the  problem  of  free  convection  from 
a  horizontal  plate  was  discussed  by  Stewartson  (32).  He  demonstrated 
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that  boundary  layer  approximations  were  no  longer  applicable.  Thus, 
purely  mathematical  treatments  of  this  problem  have  so  far  been  restricted 
to  the  one-dimensional  case.  Howard  (2)  formulated  a  one -dimensional 
problem  using  the  Boussinesq  approximations.  He  was  able  to  predict 
Townsend's  experimental  data  (3)  from  his  calculations. 

Recently,  Morton  (8)  has  discussed  the  turbulent  transport  in  a 
convection  plume  and  formulated  the  plume  problem  afresh.  His  calcula¬ 
tions  are  in  preparation. 

NUMERICAL  METHOD 

Numerical  methods  are  a  very  powerful  tool  for  solving  non-linear 
and  multi-dimensional  flow  problems.  Several  numerical  flow  calculations 
for  free  convection  problems  appear  in  the  literature.  These  calculations 
have  all  been  based  on  the  formulation  using  Boussinesq  approximations 
incorporating  the  assumption  of  constant  fluid  properties  or  constant 
eddy  dif fusivities ,  and  have  dealt  with  problems  with  well  defined 
boundaries  (within  a  rectangular  cavity,  e.g.).  In  most  procedures  a 
time  advancing  (transient)  scheme  has  been  used  to  obtain  eventually 
steady  state  solutions.  Typical  of  these  are  work  by  Wilkes  (9), 

Deardorff  (33),  and  Fromm  (34).  Torrance  (43),  in  a  recent  paper,  has 
made  detailed  comparisons  of  five  transient  type  finite-difference 
computation  schemes  (9,34,43)  for  laminar  free  convection  in  a  vertical 
axisymmetric  enclosure  with  a  small  heated  spot  centrally  located  on  the 
floor.  He  concluded  that  the  calculated  flows  for  all  methods  were 
similar  and  the  required  computer  times  were  also  of  comparable  magnitudes. 


I  i 

The  work  by  Nielsen  (13)  is  probably  the  most  directly  relevant  to 
this  thesis.  He  calculated  the  flow  field  due  to  free  convection  over  a 
heated,  horizontal,  circular  surface.  The  calculation  was  intended  for 
the  modeling  of  the  flow  fields  generated  by  a  mass  fire.  A  steady  state 
formulation  including  possible  rotation  about  the  axis  of  symmetry  was 
chosen.  The  use  of  Boussinesq  approximations  and  constant  eddy  diffusivity 
assumption  in  his  calculations  reduced  the  governing  differential  equations 
to  the  form  of  a  laminar  flow  problem  in  which  Grashof  number  is  the 
governing  parameter.  Because  of  the  rather  large  equivalent  laminar 
Grashof  number  implied  in  his  calculations  (~10^®)  through  the  choice 
of  eddy  diffusivity  (20  ft“/sec),  diameter  (4340  ft)  and  surface  tempera¬ 
ture  rise  (200°F  -  800°F) ,  the  validity  of  his  results  is  uncertain. 

His  results  showed  that  earth  rotation  has  little  effect  on  the  overall 
behavior  of  the  flow  field.  No  interpretation  or  experimental  correlation 
was  reported  in  reference  (13). 

Laminar  free  convection  due  to  a  heated  horizontal  plate  (35)  and  a 
horizontal  disk  (36)  in  a  full  space  for  small  Grashof  numbers  (<300) 
have  been  reported  in  the  recent  literature.  These  studies  were  carried 
out  using  the  Boussinesq  approximations,  consf-t..t  fiuid  properties  and  a 
steady  state  formulation.  At  the  com .us ions  of  their  report,  Kane  and 
Yang  (36)  recommended  the  study  of  cases  for  higher  Grashof  numbers  using 
a  transient  type  formulation  such  as  that  reported  in  references  (9), 

(34)  and  (43).  This  particular  work  incorporating  a  Gr^  defined 
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explicity  in  terms  of  an  equivalent  surface  temperature  rise  was  carried 
out  in  the  numerical  calculation  part  of  this  thesis.  The  details  are 
discussed  in  Chapter  VI. 


III.  DESCRIPTION  OF  EXPERIMENTAL  WORK 


The  objective  of  the  experimental  program  was  to  investigate  the 
turbulent  flow  and  temperature  fields  above  a  heated,  circular,  hori¬ 
zontal  plate  which  was  mounted  flush  with  the  ground  surface.  Both  flow 
visualization  and  flow  and  temperature  field  measurements  were  conducted 
in  this  program. 

EXPERIMENTAL  APPARATUS 

The  experimental  apparatus  for  this  study  consisted  of  three  major 
items:  the  heater  assembly,  the  traversing  mechanism,  and  the  platform 
assembly.  The  size  of  the  experimental  apparatus  was  designed  so  that 
fully  turbulent  flow  could  always  be  obtained.  Since  any  laboratory  size 
free  convection  flow  field  is  very  weak,  any  side  wind  would  easily  break 
the  flow  symmetry.  Therefore,  the  experimentation  was  conducted  indoors. 
The  overall  arrangement  of  the  apparatus  is  shown  in  Figure  la. 

Flow  visualization  is  a  useful  means  of  gaining  understanding  of  a 
flow  field.  Shadow  observation  ana  pnotography  are  well  suited  to  this 
problem.  The  required  collimacea  light  was  approximated  by  a  high 
intensity  mercury  arc  lamp  placed  at  a  large  distance  from  the  heated 
plate.  Therefore,  a  basement  laboratory  adjacent  to  a  long  hallway  was 
chosen  as  the  site  of  the  experiment.  This  location  was  also  desirable 
because  it  was  free  from  air  currents  that  woula  disturb  the  light  path. 
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Heater  Assembly 

The  beater  assembly  used  in  the  apparatus  is  shown  in  Figure  2.  A 
one-half  inch  thick  by  24  inch  diameter  aluminum  plate  was  heated  by 
radiation  from  an  electrical  heating  unit  located  directly  below  it.  The 
heater  assembly  was  mounted  in  the  center  of  an  8'  x  8'  horizontal  *  wooden 
platform  such  that  the  aluminum  plate  surface  was  flush  with  the  floor  of 
the  platform  which  simulated  the  ground  surface.  Electricity  was  supplied 
to  the  heater  unit  by  a  variable  voltage  power  transformer  (0  -  110  volt). 
By  this  means  the  power  input  to  the  heater  unit  could  easily  be  set  for 
each  experimental  run. 

Number  14  Nichrome  V  wire  was  used  for  heating.  It  was  placed  in 
equally  spaced  (0.4"  apart)  parallel  slots  cut  in  a  1/2"  thick  by  24" 
diameter  transite  board.  This  arrangement  permitted  free  lengthwise 
movement  due  to  thermal  expansion.  The  wire  was  held  in  place  by  a 
1/8"  thick  by  24"  diameter  copper  plate.  A  1/16"  thick  asbestos  sheet 
was  sandwiched  between  the  copper  and  the  transite  plates  for  electrical 
insulation.  The  heating  unit  was  placed  in  an  aluminum  pan  which  was 
insulated  on  the  inside  wall  with  a  high  temperature  insulation.  The 
aluminum  plate  was  supported  freely  by  the  heads  of  four  equally  spaced 
bolts  protruding  from  the  aluminum  pan.  The  outside  surface  of  the 
copper  plate  in  the  heating  unit  and  the  underside  surface  of  the  aluminum 

INSULAG,  k  =  0.7  But-in/hr-f t^-°F  at  1000°F  and  lower  for  lower 
temperatures.  Detailed  information  may  be  found  in  Bulletin  No. 
332-A,  Quigley  Co.,  415  Madison  Avenue,  New  York  17,  New  York. 


A  1/2"  air 
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place  were  painted  with  matte  black  high  temperature  paint, 
gap  was  maintained  between  these  two  painted  surfaces.  Due  to  the  high 
thermal  conductivity  of  copper  and  the  closely  spaced  heating  wire,  the 
temperature  of  the  copper  plate  was  very  uniform.**  This  arrangement 
permitted  indirect  heating  by  radiation  which  insured  a  uniform  heat 
flux  to  the  aluminum  plate.  The  aluminum  pan  was  supported  at  its 

jULJL 

periphery  by  the  platform.  An  insulation  jacket  was  placed  outside 
the  aluminum  pan  for  the  purpose  of  low  temperature  insulation.  The 
heat  loss  through  this  insulation  was  about  25%  of  the  total  electrical 
power  input  to  the  heater  assembly.  With  temperature  measurements  at 
various  locations  in  the  insulation  jacket,  an  energy  balance  calculation 
could  be  made  for  the  heater  assembly. 

Traversing  Mechanism 

The  traversing  mechanism  is  shown  in  Figure  la  and  lb.  It  was 
designed  to  provide  for  translation  of  a  hot-wire  probe  in  three  mutually 
perpendicular  directions  (translation  parallel  to  the  vertical  axis  of 
the  plate,  and  horizontal  translations  in  the  radial  and  the  tangential 
directions)  plus  rotation  about  the  probe  axis. 


vV 

Fuller  paint  -  flat  black  primer.  Ordinarily  used  for  painting 
fireplace  walls.  Satisfactory  for  temperatures  up  to  1000°F. 

‘'''Based  on  the  lateral  surfaca  area,  207.  of  the  heat  loss  from  the 
side  wall  of  the  insulation  jacket  was  assumed  to  be  from  the 
copper  plate.  The  radial  temperature  gradient  thus  calculated 
was  less  than  i°F/ft. 

***No.  4  vermicuiite  (a  water-repellent  masonry  fill  insulation) 
was  used  in  the  jacket  (k  =  0.55  Btu-in/hr-f t^-OF  at  200°F, 
tabulated  values  of  k  may  be  found  in  Special  Report  No,  256, 
Research  Laboratory,  Zonolite  Company,  1827  Benson  Avenue, 
Evanston,  Illinois. 
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A  traversing  h«ad  supported  by  an  8-foot  long  horizontal  bar  was  used 
co  permit  the  radial  aovaaent.  The  ends  of  this  horizontal  bar  were 
guided  by  two  vertical  supports  fastened  to  the  wooden  platform  for  the 
vertical  movement. 

The  method  of  mounting  the  hot-wire  probe  is  shown  in  Figure  3.  The 
hot-wire  probe  could  be  rotated  by  rotating  either  the  probe  holder  or 
the  probe  support.  Tangential  motion  of  the  probe  was  achieved  simply 
by  sliding  motion  of  either  the  probe  holder  or  the  probe  support. 

Small  ball  bearings  were  used  for  guiding  all  translational  movements 
except  in  the  tangential  diraction.  A  counterweight  was  used  to  balance 
the  weight  of  the  horizontal  bar  so  that  it  could  easily  be  moved 
vertically. 

Platform  Assembly 

The  platform  assembly  consisted  of  two  main  parts:  the  wooden  plat¬ 
form  and  the  screens.  The  platform  supported  both  the  heater  assembly 
and  the  traversing  mechanism.  The  floor  of  the  platform  was  made  of 
two  sheets  of  3/4"  x  4'  x  8'  plywood  placed  side  by  side.  A  circular 
opening  large  enough  for  the  heater  assembly  was  made  at  the 
center  of  the  floor.  A  small  concentric  circular  opening  was  cut  in  a 
piece  of  3/4"  x  4'  x  4'  plywood  board  bolted  below  the  floor.  The  step 
thus  formed  served  as  a  support  for  the  heater  assembly  (see  Figure  2). 
Four  pieces  of  8-foot  long  angle  iron  were  bolted  to  the  bottom  of  the 
floor  for  reinforcement.  The  floor  was  raised  to  a  level  approximately 
10  inches  above  the  floor  of  the  laboratory  and  was  firmly  supported  by 


twelve  wooden  legs . 


17 


To  minimize  the  effects  of  air  currents  inside  the  room,  the  region 
above  the  platform  was  enclosed  on  all  four  sides  by  vertical  sections  of 
ordinary  window  screen  (14  x  18  mesh  galvanized  iron)  and  by  a  horizontal 
section  of  screen  positioned  7  feet  above  the  platform.  The  screens  were 
made  such  that  the  traversing  mechanism  could  be  operated  from  outside  of 
the  screen.  Since  the  turbulent  flow  field  was  generated  from  the  hot 
place  which  was  far  from  the  screens,  and  since  the  turbulence  generated 
by  free  convection  must  greatly  exceed  that  in  a  calm  ambient  atmosphere, 
the  presence  of  the  screens  would  have  very  little  effect  on  the  turbulnece 
level  in  the  flow  field. 


FLOW  VISUALIZATION 
General  Discussion 

Due  to  the  large  size  of  the  apparatus,  many  of  the  well-established 
flow  visualization  techniques,  such  as  interferometry  or  the  Schlieren 
method,  were  ruled  out  for  this  study.  Two  basically  different  flow 
visualization  methods  were  investigated:  smoke  tracing  and  shadow¬ 
graph.  The  latter  method  was  finally  chosen  for  this  study  because  of 


A  dense  smoke  is  generated  by  pouring  FOG  JUICE  (a  petroleum 
derivative,  available  in  theatrical  supply  stores)  on  the  hot 
surface.  The  fluid  is  evaporated  upon  contact  with  the  hot  surface 
thus  forming  a  smoke  which  is  indicative  of  the  pattern  of  the 
flow  field.  The  heat  flux  available  from  the  plate  was  small 
(  —  0.25  Btu/sec-f  t2) .  Due  to  heat  removal,  the  plate  temperature 
would  drop  so  that  shortly  after  the  application  of  the  fluid,  it 
was  evaporated  in  an  erratic  manner.  This  limited  the  possible 
observation  period  to  only  a  few  minutes  at  a  time.  The  difficulty 
of  uniform  application  on  the  plate  without  disturbing  the  flow 
was  another  undesirable  feature  of  this  method. 
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the  capability  of  wide  coverage  and  unobstructed  continuous  observations. 
The  basic  requirements  for  the  shadow  method  are  a  collimated  light  beam 
and  a  projecting  screen.  A  General  Electric  B-H6  mercury  vapor  arc  lamp 
was  used  for  the  light  source.  This  was  a  quarts-capollary  enclosed  arc 
1.5  mm  x  25  mm  with  a  total  flux  of  60,000  lumens.  The  lamp  dissipated 
900  watts  and  was  cooled  by  a  30  psig  air  supply.  In  order  to  approximate 
a  collimated  light  beam,  this  light  source  was  placed  at  the  far  end  of 
a  hallway  in  the  basement  laboratory  («s#165  feet  from  the  hot  plate).  A 
translucent  viewing  screen  made  of  a  sheet  of  ordinary  tracing  vellum 
was  placed  normal  to  the  light  path  about  20  feet  on  the  other  side  of 
the  hot  plate.  Due  to  the  light  transmission  of  the  screen,  it  was  found 
quite  satisfactory  to  view  the  shadow  from  hehind  the  screen. 

For  photography,  a  piece  of  household  curtain  material  (aluminum 
coated  lining)  was  hung  on  the  wall  normal  to  the  light  path  about  27 
feet  on  the  other  side  of  the  hot  plate.  Photographs  were  taken  in 
front  of  the  screen  at  an  angle  just  large  enough  not  to  block  the  light 
path  (Negative  Tri  X  Pan  35  mm  film  was  exposed  f2  at  1/30  second  and 
developed  in  Acufine  for  an  ASA  rating  of  1200). 

Flow  Field  Description 

The  flow  field  about  the  heated  plate  was  viewed  using  the  shadow 
method.  Figure  A  shows  the  shadow  pictures  taken  at  two  plate  mean 
temperatures  by  photographing  the  reflection  of  the  shadows  from  the 
aluminum  coated  lining  screen.  These  pictures  show  that  the  overall 
flow  field  consists  of  two  major  regions  which  are  separated  by  a  fairly 
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well  defined  boundary.  Above  Che  plate  and  inside  this  boundary  is  a 
region  here  called  the  primary  flow  region  in  which  there  is  a  dominant 
upward  moving  stream  of  eddies.  The  upward  flowing  column  is  fully 
turbulent  and  is  constricted  about  the  axis  of  symmetry  (the  centerline) 
of  the  plate.  It  may  be  seen  from  Figure  4  that  the  flow  boundary  is 
more  sharply  defined  for  higher  plate  temperature  (or  equivalently  higher 
heat  flux).  This  phenomenon  has  been  repeatedly  observed. 

Outside  the  boundary  is  the  region  of  secondary  flow  which  moves 
across  the  boundary  into  the  primary  flow  region  to  replace  the  upward 
moving  hot  air.  The  temperature  in  the  secondary  flow  is  substantially 
uniform.  Because  of  this,  it  is  impossible  to  visualize  in  the  shadow 
picture.  The  shape  of  the  flow  region  boundary  together  with  knowledge 
of  the  velocities  from  direct  observation  suggest  that  the  secondary 
flow  (the  indraft)  velocity  is  largest  within  a  few  inches  above  the 
ground. 

It  was  observed  that  the  flow  field  could  be  upset  by  even  a  small 
ambient  disturbance.  This  is  to  be  expected,  in  view  of  the  very  weak 
velocity  field.  Careful  observation  of  the  shadow  revealed  an  interesting 
feature  of  the  disturbed  flow  field.  In  a  closed  room  of  still  air,  the 
flow  field  remains  symmetric.  If  a  small  disturbance  is  introduced, 
e.g.,  by  placing  a  small  piece  of  curled  paper  near  the  edge  of  the  plate, 
the  flow  field  sometimes  is  disturbed.  A  clearly  defined  continuous 
whirling  column  of  finite  width  may  be  seen,  as  illustrated  in  Figure  5. 
This  whirl  starts  from  the  plate  surface  and  eventually  penetrates  upward 


for  several  heater  place  radii.  In  shape  chis  whirl  say  be  scraighc  or 
curved,  wide  or  narrow.  Ic  may  stand  still  or  wander.  Ics  life  may 
vary  from  a  fraction  of  a  second  to  a  few  seconds.  For  large  disturbances 
such  as  bulk  motion  of  the  room  air,  che  occurrence  of  this  whirl 
becomes  much  more  frequent  and  the  whirl  width  may  be  as  wide  as  2  inches. 
Further  studies  of  this  phenomenon  will  certainly  be  of  interest. 

For  the  purpose  of  this  experimental  study,  flow  field  symmetry  is 
essential.  In  order  to  minimize  any  room  disturbances  which  might 
break  the  flow  symmetry,  the  screens  described  in  the  previous  section 
were  installed.  Shadow  observations  showed  that  these  screens  did 
indeed  stabilize  the  flow  field. 

1  le  turbulence  in  the  primary  flow  region  is  of  course  another 
interesting  feature  of  the  flow  field.  It  is  more  appropriate  to  describe 
the  turuulence  in  conjunction  with  the  hoc-wire  data  1 5  the  discussion 
of  the  turbulence  will  be  postponed  until  Chapter  V. 

HEATER  PERFORMANCE 

In  free  convection  flows,  che  piate  temperature  and  the  heat  transfer 
are  related.  The  performance  of  the  heater  can  be  characterized  by 
either  the  mean  temperature  difference  between  the  plate  and  ics  ambient 
or  the  mean  heat  flux. 

Plate  Temperature 

The  surface  temperature  distribution  of  the  heated  aluminum  plate 
was  measured  by  thirteen  thermocouples  installed  in  the  four  radially 


milled  sloes  on  che  underside  of  the  place  (Figure  6).  These  therno- 
ccuples  were  cade  from  Kucher  24  "high  accuracy"  copper-con* cancan  wire 
(L  £  H  Wire  No.  24-55-11,  error  1/4°F  ac  175°?  with  reference  junction 
at  32°F),  and  were  connected  to  a  24  channel  switching  box  which  was 
directly  connected  to  a  potentiometer.  Xly  time  mean  temperatures  were 
oeasured.  Figure  7  shows  the  measured  plate  surface  temperature  distri¬ 
butions  for  two  mean  plate  temperatures.  It  is  noticed  that  the  plate 
temperature  is  quite  uniform  in  the  central  ragion  and  drops  off  slightly 
toward  the  edge  of  the  plate.  This  fact  is  expected,  and  attributed  to 
the  more  effective  cooling  by  the  cold  air  near  the  edge  of  the  plate 
(considering  the  fact  of  uniform  heat  flux  input  to  the  plate).  Assuming 
chat  the  heat  flux  is  proportional  to  the  4/3  power  of  the  plate  tempera¬ 
ture  rise  above  the  ambient  (27),  this  temperature  non-uniformity  repre¬ 
sents  a  maximum  of  3%  local  heat  flux  variation.  As  far  as  concerns  the 
boundary  condition  for  the  overall  free  convection  flow  field,  this 
slight  temperature  non-uniformity  is  insignificant. 

Energy  3alance 

Twenty  No.  24  iron-cons  cancan  thermocouples  (L  &  N  Wire  No.  24-50-39) 
were  installed  in  che  insulation  jacket  of  heater  assembly  for  energy 
balance  measurements  (Installation  details  are  given  in  Figure  2).  These 
thermocouples  were  also  connected  to  a  24  channel  switching  box  through 
which  any  one  of  them  could  be  connected  to  a  potentiometer  for  reading. 

The  total  energy  loss  from  the  heater  assembly  is  the  sum  of  energy 
lost  to  the  surrounding  air  by  radiation  from  the  aluminum  plate  (Qp), 


by  convection  and  radiation  from  the  asbestos  ring  (Qa)»  by  convection 
and  radiation  from  the  floor  of  the  platform  (Q^),  and  by  conduction 
through  the  insulation  jacket  (Q^). 

The  net  convective  heat  transfer  from  the  heated  plate. to  the  air 
above  (Qc)  is  computed  as  follows: 

Qc  '  «1„  -  9i3ss  -  3-413  El  -  <Qp  +  Q,  1-  Qf  +  Qj) 

where  E  *  voltage  drop  across  the  heater,  volt 

I  *  current  passing  through  the  heater,  ampere 
Q  *  heat  transfer  rate,  Btu/hr 

Energy  balances  were  made  for  two  heater  power  settings.  The  detailed 
calculations  of  these  terms  are  given  in  Appendix  B.  A  summary  of  all 
the  loss  terms  is  tabulated  below: 


Table  1.  Energy  balance  terms. 


Run 

No. 

E 

volt 

H 

°F 

Op 

Qin 

3/hr 

B/Sr 

Qa 

B/hr 

Qf 

B/hr 

B/ir 

Qc 

B/hr 

1 

132 

8.6 

520 

77 

3880 

229 

525 

241 

347 

2538 

2 

108 

410 

77 

2615 

133 

359 

114 

240 

1769 

It  is  clear  that  the  loss  terras  are  all  small  relative  to  Qc.  Therefore, 
the  percent  accuracy  of  Qc  is  better  than  any  of  the  heat  losses.  As  a 
check,  a  direct  energy  balance  measurement  was  made  by  blocking  off  the 
heat  dissipation  from  the  aluminum  plate.  In  this  manner,  the  energy 
input  to  the  heater  assembly  at  any  given  plate  temperature  would 
represent  the  net  heat  loss  from  the  assembly  (less  the  radiative  heat 
loss  from  the  aluminum  plate).  The  results  of  this  direct  measurement 
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and  those  calculated  based  on  temperature  measurements  are  shown. in 
Figure  8. 

Heat  Transfer  Correlation 

An  examination  of  the  significant  governing  variables  for  the  heat 
transfer  phenomenon  leads  to  the  following  relation: 

Nu  -  f  (Pr ,  Gr „  0Af) 
where  Nu  =  hD/k,  Nusselt  number, 

Pr  ■  **Cp/k,  Prandtl  number, 

Gr  =  g£A0B3  p2/*2,  Grashof  number, 
and  h,  D,  k,  p,  g,  &  AP,  and  p  are,  respectively,  the  average  heat 
transfer  coefficient,  plate  diameter,  thermal  conductivity,  dynamic 
viscosity,  gravitational  acceleration,  coefficient  of  volumetric  expansion, 
temperature  difference  between  the  plate  mean  and  the  ambient,  and  density. 
For  weakly  buoyant  flows  (4,27,37,38)  such  as  in  free  convection,  the 
parameter  0A0  is  not  significant  enough  to  affect  the  flow  by  itself  and 
is  absorbed  in  Gr.  This  can  also  be  seen  from  the  governing  differential 
equations  for  small  density  variations  (see  Chapter  VI). 

For  free  convection,  the  flow  becomes  fully  turbulent  when  the 
Grashof  number  is  very  large  (>10^).  This  can  be  achieved  most  effectively 
by  making  D  sufficiently  large.  Since  the  heat  flux  of  a  very  large 
surface  can  only  be  finite  and  cannot  be  affected  by  changing  its  size, 
we  assume  that  the  heat  flux  is  independent  of  the  size  D,  i.e., 

Nu  «  (Gr)1/3 


Furthermore,  experimental  correlations  (26,27)  indicate  that  the  primary 
effect  of  Pr*  can  be  taken  into  account  by  grouping  it  with  Gr  to  form  a 
Rayleigh  number  Ra,  i.e., 

Nu  «  constant  x  Ra^^ 

where  Ra  *  PrGr.  Using  the  property  values  of  air  at  the  ambient  tempera- 

2 

ture,  and  the  definition  of  heat  flux  (q,  Btu/hr-ft  ) 

q  =  hA# 

Nu  and  Ra  for  the  two  runs  in  Table  1  may  be  calculated.  The  relevant 
values  are  given  below: 


Table  2.  Dimensionless  heat  transfer  parameters. 


Run 

No. 

q,  Btu/ 
hr-ft2 

h,  Btu/ 
hr-ft  -°F 

Pr 

Nu 

Gr  x  10" 10 

Ra  x  10'10 

800 

1.81 

0.72 

244 

0.772 

0.556 

LiJ 

560 

1.68 

0.72 

226 

0.573 

0.413 

The  above  data  are  plotted  in  Figure  9  and  are  found  to  be  in  good 
agreement  with  experimental  results  of  Fishenden  and  Saunders  (26).  One 
may  observe  that  the  result  Nu  oc  Ra^3  is  consistent  with  the  present 
data,  as  expected.  The  proportionality  constant  was  calculated  to  be 
0.141,  so  that: 

Nu  =  0.141  x  Ra1/3.  (3-1) 


The  Prandtl  number  for  air  is  substantially  a  constant  in 
the  range  of  temperature  in  this  investigation. 
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VELOCITY  MEASUREMENT 
General  Discussion 

We  are  concerne  *.  with  the  velocity  measurement  in  the  turbulent 
free  convection  regime.  In  a  flow  field  of  this  type  local  mean  velocity 
varies  up  to  3  ft/sec,  but  the  velocities  over  most  of  the  flow  field  are 
only  on  the  order  of  a  fraction  of  a  foot  per- second.  The  turbulence  in 
the  flow  field  is  characterized  by  low  frequency,  high  amplitude  fluctua- 
tions  (see  Chapter  V  for  detailed  discussions).  Very  little  has  been 
reported  in  the  literature  concerning  this  type  of  velocity  measurement. 
Tritton  (28,39)  used  a  quartz  fiber  to  measure  the  mean  velocities  in  the 
boundary  layer  region  above  an  almost  horizontal  heated  plate.  He  used 
a  telescope  to  observe  the  deflection  of  the  cantilevered  quartz  fiber 
in  the  flow  field.  Since  there  is  no  easily  recorded  signal,  this  method 
is  not  considered  suitable  for  large  number  of  measurements.  Very 
recently,  the  laser  doppler  technique  (14,15)  has  been  successfully 
applied  to  turbulent  flow  field  measurement.  Aside  from  the  obvious 
advantage  of  undisturbed  measurement,  it  is  possible  to  use  this  technique 
to  measure  both  the  instantaneous  and  mean  velocity  component  in  a  desired 
direction.  No  attempt  was  made  to  use  this  technique  for  the  present 
studies  since  the  measurements  for  this  experiment  were  well  underway 
before  the  possibilities  of  using  such  a  technique  became  apparent. 

Tne  hot-wire  technique  is  very  useful  for  turbulence  measurements. 
However,  in  slow  velocity  measurements,  it  suffers  severely  from 
orientation  effects  due  to  the  free  convection  field  generated  by  the 


26 


wire  itself.  Nevertheless,  the  ease  of  meeting  the  time  response 
requirement,  ease  of  recording,  and  the  immediate  availability  of  the 
required  electronic  components  prompted  the  decision  to  use  the  hot-wire 
method  for  this  experiment.  In  order  to  minimize  the  effects  due  to 
heat  generated  in  the  wire,  earlier  attempts  were  centered  at  making 
the  sensing  wire  as  small  as  possible.  The  smallest  wire  used  was 
0.00012"  diameter  Pt-107»  Rh  wire.  Unfortunately,  wire  of  this  size  did 
not  survive  the  vibration  resulting  from  large  probe  traveling.  Moreover, 
as  long  as  the  wire  is  at  a  higher  temperature  than  the  ambient,  free 
convection  effects  are  always  present. 

Collis  (40,41)  in  his  reports  on  forced  convection  of  heat  from 
cylinders  at  low  Reynolds  numbers  (Re  =  dV/»/*)>  gave  the  following 
criterion: 

1/3 

Free  convection  effects  are  negligible  when:  Re  >  Gr 

This  expression  implies  that  free  convection  effects  are  unimportant  when 
the  local  flow  velocity  undisturbed  by  the  wire  is  larger  than  (g£Atfv)^. 
It  shows  also  that  free  convection  effects  are  independent  of  the  wire 
diameter  and  that  they  can  be  reduced  by  operating  the  wire  at  small 
values  of  Ad,  the  temperature  difference  between  the  wire  and  the  ambient. 

The  resistance  and,  hence,  the  calibration  of  a  thin  wire  tends  to 
drift  due  to  self-straining.  This  situation  may  be  improved  by  choosing 
a  large  wire  diameter  and  by  properly  soldering  the  wire  to  its  supports. 
Large  wire  diameter  is  desirable  for  its  strength  but  it  reduces  the  wire 
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resistance  and  hence  the  output  signal  (the  voltage  drop  across  the  wire) 
for  fixed  temperature  operation.  As  a  reasonable  compromise,  the  hot¬ 
wire  chosen  was  bare  platinum  wire  0.001"  diameter  by  <5.'5V  long.  This 
kind  of  wire  was  found  to  be  quite  satisfactory  in  terms  of  strength, 
repeatability  of  results,  sensitivity,  directional  sensitivity,  and 
fast  time  response  (see  Chapter  V  for  more  discussions  on  the  choice  of 
this  wire  length).  The  time  constant  for  this  wire  was  calculated  to  be 
less  than  one  millisecond  for  constant  temperature  operations  (42).  In 
order  to  minimize  the  effects  due  to  the  free  convection  field  induced 
by  the  hot-wire,  a  low  wire  operating  temperature  (i.e.,  a  small  overheat 
ratio)  was  used. 

Hot-Wire  Calibration 

A  TSI  model  1010  constant  temperature  anemometer  was  used  for 
generating  the  hot-wire  signals.  In  order  to  account  for  the  free  convec¬ 
tion  effects,  extensive  calibrations  were  made  in  which  both  the  wire 
orientation  and  the  velocity  vector  were  taken  into  consideration. 

A  calibration  tunn^i.  was  designed  to  meet  the  requirements  of 
variable  flow  velocity,  orientation,  and  air  temperature.  Figure  10 
shows  a  schematic  diagram  of  this  tunnel.  The  tunnel  test  section  was 
an  enclosed  4"  T.D,  alumin -n  can,  the  axis  of  which  could  be  rotated 
through  an  angle  of  180  degrees  on  a  supporting  bracket.  An  air  jet 
issued  from  a  1"  diameter  orifice  with  a  rounded  inlet.  Six  sheets  of 

* 

Thermo-Systems,  Inc.,  2500  Cleveland  Avenue  North,  St. 

Paul,  Minnesota  55113. 


28 


spaced  100  x  100  mesh  bronze  wire  strainer  cloth  were  placed  upstream  of 
the  orifice  to  straighten  the  flow.  The  hot-wire  probe  was  inserted  into 
the  tunnel  through  a  removable  support  plug,  when  the  wire  was  in  place 
it  was  about  1/2"  downstream  of  the  orifice  exit  plane.  The  probe  is 
shown  approximately  to  scale  in  Figure  3. 

The  uniformity  of  the  velocity  distribution  in  the  jet  was  checked 
by  traversing  a  0.2"  long  hot-wire  probe  along  the  jet  diameter.  The 
velocity  distribution  in  the  center  region  of  the  jet  was  essentially 
uniform.  \  drop  of  approximately  1%  of  the  hot-wire  bridge  voltage  output 
was  observed  at  957„  radius.  With  the  hot-wire  placed  in  the  center  of 
the  jet  during  calibration,  the  total  length  of  the  hot-wire  was  exposed 
to  a  uniform  velocity  stream. 

The  velocity  at  the  orifice  could  be  varied  from  0  to  3  ft/sec  by  a 
control  valve  which  was  located  between  the  outlet  of  a  pressure  regulator 
and  the  inlet  of  a  Manostat  flow  meter.*  This  meter  was  used  for 
measuring  velocities  below  2  ft/sec.  A  water  U-tube  manometer  was  used 
for  measuring  higher  velocities.  Figure  11  shows  the  calibration  curves 
for  the  flow  meter  and  the  U-tube  expressed  in  terms  of  the  orifice  exit 
velocity  at  a  standard  condition  (70°F  and  30"  Hga)  versus  the  meter  or 
U-tube  readings.  The  temperature  of  the  jet  could  be  varied  through  the 
use  of  a  100  watt  cartridge  type  electric  heater  which  could  heat  the  air 
up  to  160°F  at  an  orifice  exit  velocity  of  2  ft/sec. 

«Kp 

Manufactured  by  Manostat  Corporation,  Mew  York. 
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The  desired  wire  operating  temperature  was  determined  by  comparing 
the  calibration  at  several  wire  temperatures.  It  was  found  that  operating 
the  wire  at  approximately  300°F  was  quite  satisfactory  for  the  present 
applications.  At  this  wire  temperature,  a  reasonably  good  wire  sensitivity 
could  be  obtained  and  the  free  convection  effect  was  found  negligible 
for  velocities  greater  than  0.15  ft/sec. 

Figure  12  shows  a  set  of  hot-wire  calibration  curves.  This  set  of 
curves  consists  of  the  calibrations  of  a  horizontal  wire  at  constant 
wire  operating  temperature  (~300°F  or  4.50  fi  wire  operating  resistance) 
but  at  three  different  ambient  air  temperatures.  The  calibration  data 
were  fitted  into  the  following  functional  form  by  the  method  of  least 
squares : 

E2  =  A  +  BW1/2  (3-2) 

where  E  is  the  hot-wire  bridge  voltage  output,  W  the  velocity,  and  A  and 
B  are  constants  depending  on  the  ambient  temperature  through  their 
dependence  on  fluid  properties.  Direct  comparison  of  the  data  for  air 
blowing  upward  and  blowing  downward  shows  that  the  free  convection  effect 
is  indeed  small  for  W<0.15  ft/sec. 

For  small  velocities  on  the  order  of  0.15  ft/sec,  another  set  of 
calibration  curves  for  the  hot-wire  was  prepared  (Figure  13).  It  is 
seen  that  the  velocity  vector  of  a  laminar  flow  can  be  determined  once  a 
set  of  two  hot-wire  outputs  (for  wire  horizontal  and  vertical)  at  the 
same  location  is  obtained.  This  set  of  calibration  curves  was  used  only 
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for  reducing  the  indraft  velocity  data  for  which  the  flow  was  essentially 
laminar  and  hence  no  non-linear  averaging  process  was  involved  as  would 
be  the  situation  for  turbulent  flows. 

It  is  remarked  that  for  small  velocity  and  laminar  flow  such  as  the 
indraft  velocity  measurements,  wire  orientation  effects  are  important 
and  therefore  were  considered  in  determining  the  velocity  vector.  For 
the  turbulent  flow  regime  of  this  study  however,  the  wire  output  is 
substantially  independent  of  orientation  because  of  the  relatively  high 
fluctuating  velocities. 

Further  details  on  hot-wire  data  taking  and  the  response  equation 
will  be  discussed  where  appropriate  in  Chapter  V. 


IV.  TEMPERATURE  FIELD  AND  ANALYSIS 


The  overall  temperature  field  may  be  divided  into  two  main  regions: 
the  region  very  near  the  plate  (within  1/2"  of  the  plate"  and  the  region 
above.  In  the  region  near  the  plate,  a  large  temperature  gradient 
exists.  The  major  portion  of  the  overall  temperature  drop  takes  place 
in  this  region.  In  the  region  above,  both  the  air  temperature  and  the 
temperature  gradient  are  much  smaller  than  those  in  the  near-plate-region. 
For  this  reason,  the  measurements  for  these  regions  were  taken  at  separate 
times  with  different  thermocouple  mounting  methods.  They  will  be  discussed 
under  separate  headings. 

Air  temperatures  in  both  of  these  regions  were  measured  with  a  No.  34 
gage  copper-constantan  thermocouple.  This  thermocouple  was  calibrated 
against  a  pair  of  certified  Pt  and  Pt-Rd  thermocouples  in  an  oil  bath 
(see  calibration  curve  in  Figure  14).  The  thermocouple  data  output  was 
displayed  on  a  BRUSH  continuous  chart  recorder  via  the  amplifier  output 
circuit  of  a  Hewlett  Packard  Type  413A  DC  null  voltmeter.  The  time 
average  of  the  thermocouple  output  (time  constant  ~ 1  sec)  was  obtained 
by  use  of  a  large  time  constant  (~1  minute)  integrating  circuit  similar 
to  those  shown  in  Figure  15.  The  temperature  at  each  measured  location 
was  recorded  on  the  chart  for  a  period  of  two  to  five  minutes  and  the 
steady  state  thermocouple  output  was  used. 


MEAN  TEMPERATURE  FIELD  ABOVE  TK£  PLATS 

The  mean  temperature  field  was  o&isured  simultaneously  with  the 
ve’ocity  field.  The  thermocouple  was  Bounced  on  che  hoc-wire  probe 
(Figure  16)  and  traversed  with  it.  Since  cemperacure  is  a  scalar  quantity 
and  the  response  of  a  thermocouple  is  nearly  linear  with  temperature  in 
the  range  cf  interest,  che  temperature  field  measurement  is  a  relativ-iy 
simple  task  compared  with  that  uf  the  velocity  field.  Several  runs  were 
made  for  each  of  the  two  mean  plate  temperatures  investigated. 

In  order  to  determine  the  symmetry  of  the  temperature  field,  measure¬ 
ments  were  made  at  three  locations  (on  the  same  radius  but  spaced  at  90° 
apart).  Data  showed  that  discrepancies  among  "nese  three  readings  were 
generally  very  small.  Tne  largest  (observed  only  occasionally)  were 
less  than  3%  of  he  average  reading.  The  discrepancy  was  believed  to  be 
mainly  due  to  the  large  period  .=?riacion  of  the  air  temperature  in  the 
laboratory  room  and  to  a  less  extent  due  to  plume  leaning.  Shadow  viewir. 
showed  that  the  plume  was  almost  always  straight  up  even  when  che  labora¬ 
tory  door  was  wide  open,  and  that  the  ouration  of  occasional  plume  leaning 
was  too  short  to  affect  che  mean  temperature  readings  significantly.  In 
view  of  the  data  scatter  ( >  ±  3°F)  shown  in  Figures  17,  18  and  19, 
asymmetry  of  the  data  was  considered  to  be  relatively  insignificant. 

For  each  plate  temperature  the  centerline  temperature  profile  and 
three  horizontal  temperature  profiles  were  obtained  (Y  =  3",  6",  12" 
abov  the  plate).  Figure  17  shows  che  centerline  temperature  profiles 
for  both  plate  temperatures  investigated.  Figures  18  and  19  show  the 


horizontal  temperature  profiles  measured  for  the  two  plate  temperatures. 
In  these  figures,  each  point  represents  the  average  of  six  to  ten  mean 
temperature  readings.  A  line  is  drawn  through  each  average  temperature 
to  indicate  the  error  of  *  one  standard  deviation.  It  is  interesting  to 
note  that  the  temperature  drops  off  almost  linearly  vith  radius  ia  the 
region  near  the  centerline. 

Figures  18  and  19  show  that  lover  average  air  temperature  and 
relatively  flat  horizontal  temperature  profiles  are  associated  with  lower 
plate  temperature.  These  fiat  average  temperature-  profiles  imply  a 
larger  plume  width  in  agreement  with  Mortor/s  analysis  (38)  that  the 
width  of  a  weakly  buoyant  plume  is  larger  than  that  of  a  strongly  buoyant 
plume.  Because  of  the  limited  amount  of  data,  quantitative  determination 
of  plume  width  was  not  practical. 

MEAN  TEMPERATURE  NEAR  THE  PLATE 

The  temperatures  near  the  plate  surface  were  measured  with  the  same 
thermocouple  (34  gauge  Cu-conscantan  wire)  as  was  used  for  the  tempera¬ 
ture  field  probing;  however,  a  different  mounting  method  was  employed 
(Figure  20).  Temperatures  were  measured  at  three  radial  locations  (X  = 
0",  6",  12")  for  each  vertical  position  Y.  Because  of  the  large  tempera¬ 
ture  gradient  existing  near  the  piate,  a  small  increment  in  elevation  Y 
was  necessary.  This  small  increment  was  obtained  through  the  use  of  a 
micrometer  attached  to  one  end  of  the  horizontal  bar  of  the  traversing 
mechanism.  By  advancing  the  micrometer  screw  againsc  a  hard  stationary 


surface,  the  whole  traversing  mechanism  could  be  raised  to  any  desired 
height  with  respect  to  the  plate  in  increments  of  0.001“. 

Figures  21  and  22  show  the  temperature  distribution  in  the  region 
near  the  plate  surface  for  the  two  place  temperatures  investigated 
(ip  *  520°F  and  410°F).  these  temperatures  were  obtained  by  correcting 
the  raw  data  for  losses  due  to  radiation  and  conduction  (see  Appendix  C 
for  details  of  calculations).  One  may  observe  two  interesting  facts  from 
these  temperature  data.  The  major  portion  of  the  overall  temperature 
drop  takes  place  within  approximately  0.1"  of  the  plate.  This  region  is 
thicker  for  the  higher  place  temperature. 

The  surface  temperature  gradient  is  larger  for  the  higher  plate 
temperature  and,  because  of  the  more  effective  cooling  near  the  plate 
edge,  this  temperature  gradient  increases  with  the  radius.  This  cooling 
effect  is  also  reflected  in  the  plate  radial  temperature  distribution  as 
shown  in  Figure  7. 

The  temperatures  within  0.1"  of  the  plate  were  calculated  based  on 
the  average  heat  flux  and  molecular  conduction  only  (solid  lines  in 
Figures  21  and  22).  These  calculations  show  that,  very  close  to  the 
plate  surface,  the  measured  temperatures  are  nearly  the  same  as  those 
based  on  conduction  calculations.  This  region  is  within  approximately 
0.02"  of  the  plate  surface.  We  will  follow  Townsend  (3)  and  call  this 
region  the  conduction  layer. 
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THERMAL  DIFFUSIVITY 

Some  estimates  on  the  turbulent  transport  coefficient  based  on  the 
available  temperature  data  are  made  here.  Adopting  the  thermal  diffusivity 
concept,  the  steady  state  energy  equation  may  be  written  as  follows: 


V  -VB  =  (  v29  +  Ve  -V# 

n  n 


or. 


H  ‘  - - - 


& 


where  V,  B,  *  ,  and  v  are,  respectively,  the  velocity  vector,  mean 

n 

temperature,  thermal  diffusivity  and  gradient  operator.  is  calculated 
using  the  following  procedure:  We  proceed  arbitrarily  to  neglect  the 
spatial  variation  of  .  This  gives: 


*H  V"8 


(4-1) 


The  value  of  «  calculated  from  Equation  (4-1)  are  then  used  to  estimate 
the  magnitude  of  the  neglected  term  Vt-^'VO  to  see  if  the  procedure  is 
justified . 

Based  on  the  mean  temperature  data  (Figure  18)  the  values  of  VO 

2  *♦ 

and  ve  were  obtained  graphically.  Using  the  data  of  V  from  Chapter  V, 

the  values  of  were  calculated  from  Equation  (4-1)  for  the  case  8^  = 
520°F.  The  values  of  and  for  several  locations  in  the  flow  field 
are  given  in  Table  3. 

The  estimates  in  Table  3  show  that  the  values  of  are  nearly 
constant  in  the  major  portion  of  the  primary  flow  region.  is  large 
near  the  plate  and  is  relatively  small  at  the  location  (6,3)  which  is 
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near  cne  boundary  of  the  primary  flow  region.  The  magnitudes  of  the 
neglected  term  are  estimated  and  are  also  give  in  Table  3.  It  is 

seen  cnat  they  are  much  smaller  than  the  typical  values  of  the  term 
Vd .  Neglecting  of  makes  the  values  of  somewhat  coo  high. 

Table  3.  Thermal  diffusivity  estimates  (*D  *  520°P). 


Location 

in. 


Velocity 
f  ps 


.93  162.5 

1.35  152 

0  1.90  140 

2.22  125.5 


.78  137 

1.06  120 


.41  .45  115.3 
.26  .77  98.7 


‘H 

°F/sec 

°F/sec 

f t^/sec 

93 

11 

.066 

88 

5 

.053 

68 

.7 

.04 

52 

.1 

.037 

26.2 

2.8 

.052 

28.4 

1.4 

.043 

4.7 

.1 

.01 

5.3 

.3 

.037 

V.  RESULTS  AND  ANALYSIS  0?  HOT-WIRE  DATA 


INTRODUCTION 

The  experimental  measurements  of  the  flow  above  a  heated  horizontal 
circular  plate  are  described  here.  The  mean  velocities  and  turbulence 
kinetic  energies  are  of  primary  interest.  Using  the  turbulence  measure¬ 
ments,  the  eddy  diffusivity  will  be  estimated. 

As  will  be  shown,  the  mean  flow  velocities  are  small  ('-1  ft/sec)  and 
the  turbulence  velocities  are  on  the  same  order  of  magnitude  as  the  mean. 
Due  to  the  large  fluctuations,  coaveatiocal  experimental  methods  and 
attendant  mathematical  simplifications  for  turbulence  measurements  fail. 

No  data  directly  relevant  to  tne  present  investigation  have  been  published. 

In  any  particular  geometric  space,  the  desired  flow  properties  are 
all  determined  by  the  probability  density  function  over  four-dimensional 
velocity-temperature  space.  The  data  reduction  procedure  is  based  on 
the  idea  that  the  essential  features  of  this  probability  density  function 
are  determined  by  the  one-point  probability  density  distribution  of  hot¬ 
wire  output  voltage  with  wire  orientation  as  a  parameter.  This  informa¬ 
tion  is  approximated  by  a  set  of  averages  of  hot-wire  output  repeated  at 
several  orientations. 

A  constant  temperature  hot-wire  anemometer  was  used  for  this  experi¬ 
ment.  The  one-point  averaged  values  of  the  hot-wire  output  and  selected 
low-order  moments  were  obtained  at  selected  wire  orientations.  The 


quantities  of  interest  were  then  inferred  from  knowledge  of  these  moments 


ana  a  sec  of  numerically  simulated  hoc-wire  data.  In  the  following 
section,  the  hot-wire  data  from  direct  experimental  measurements  will 
be  suosaarized  acd  a  numerical  simulation  of  the  turbulence  based  on  an 
idealized  computation  model  will  be  described.  A  data  reduction  proce¬ 
dure  using  the  above  results  is  then  discussed  and  the  reduced  data 
presented.  Finally,  the  justification  of  the  numerical  model  and  the 
assumptions  incorporated  in  it  will  be  given. 

HOT-WIRE  DATA  SUMMARY 
Ins  trumenta  t ion 

The  complete  instrumentation  setup  is  shown  in  Figure  15.  It 
consists  primarily  of  a  constant  temperature  hot-wire  anemometer,  two 
linear izcrs,  the  integrating  circuits  and  the  recording  devices.  One- 
point  time  averages  of  the  quantities  shown  in  Figure  15  (i.e.,  E,  E, 

E^,  etc.)  are  the  basic  data  to  be  taken. 

A  0.001"  diameter  x  0.5"  long  bare  platinum  wire  was  chosen  as  the 
sensing  element.  Shadow  observations  showed  that  the  eddy  size  was  on 
the  order  of  1",  therefore,  a  1/2"  long  wire  would  be  a  reasonable  choice. 
In  order  to  determine  the  length  effect  on  the  resolution  of  turbulence 
sensing,  a  0.2"  long  wire  was  also  made.  These  two  wires  were  operated 
at  the  same  wire  temperature  (approximately  300°F)  and  were  placed  in 
the  flow  field  6"  above  the  plate.  The  hot-wire  outputs  were  recorded 
on  an  AMPEX  Model  FR-1200  tape  recorder  at  a  recording  speed  of  1-7/8 
inches  per  second  (ips).  Comparisons  were  made  by  playing  the  tapes 
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back  at  60  ips  and  analyzed  with  a  Hewlett  Packard  Model  302A  wave 
analyzer.  As  shown  in  Figure  23,  there  was  no  significant  difference 
between  the  two  wire  outputs  in  either  frequency  or  amplitude  of  the 
fluctuations.  It  was  therefore  concluded  that  for  this  investigation 
wire  length  within  the  range  studied  was  not  an  important  parameter. 

The  longer  wire  (0 . 5‘*  long)  was  used  for  all  measurements  because 
of  the  desired  larger  voltage  output  or  better  sensitivity. 

Characteristic  Quantities  of  the  Turbulence 

Figure  24  shows  some  typical  hot-wire  output  traces  for  both 
horizontal  and  vertical  orientation  of  the  wire  at  various  distances 
above  the  center  of  the  plate.  Figure  25  shows  hot-wire  output  traces 
for  a  horizontal  wire  at  five  radial  locations  6"  above  the  plate.  This 
figure  shows  that  the  turbulence  degenerates  with  increasing  radius. 

At  a  radius  of  9"  the  traces  indicate  intermittent  turbulence.  At  6" 
radius,  the  large  amplitude  and  the  erratic  excursion  of  the  traces  suggest 
that  the  hot-wire  is  located  approximately  on  the  boundary  of  the  primary 
flow  region.  Both  Figures  24  and  25  show  that  the  amplitude  of  the 
fluctuations  are  on  the  order  of  one  volt.  This  corresponds  to  velocity 
fluctuations  on  the  order  of  1  ft/sec  (Figure  13),  which  is  of  the  same 
order  of  magnitude  as  the  mean  flow  velocity. 

Since  the  zero  velocity  position  of  these  traces  depends  on  the 
local  temperatures  which  vary  from  location  to  location,  it  is 
impractical  to  mark  these  traces  with  the  corresponding  velocities. 
Therefore,  by  themselves  these  hot-wire  traces  (Figures  24  and  25) 
demonstrate  only  the  qualitative  nature  of  the  turbulence. 
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The  use  of  a  wave  analyzer  was  originally  intended  to  obtain  the 
characteristic  frequency  of  the  turbulence.  Analysis  of  the  data  showed 
that  the  main  contribution  of  the  frequency  was  in  the  range  below  10  Hz. 
At  frequencies  below  25  Hz  the  amplitude  dropped  off  rapidly  to  less 
than  one  percent  of  the  value  at  1  Hz  and  became  negligibly  small  at 
higher  frequencies  (e.g.,  Figure  23).  Due  to  the  poor  low  frequency 
capability,  the  wave  analyzer  could  not  be  used  to  obtain  a  meaningful 
time  scale  for  the  turbulence.  Therefore,  the  autocorrelation  function 
of  the  recorded  turbulence  data  was  obtained  using  a  PAR  Model  101 
correlation  function  computer  (manufactured  by  Princeton  Applied  Research 
Corporation) . 

The  autocorrelation  function  of  the  hot-wire  bridge  voltage  output 
E  is  defined  as: 

*(t)  =  lim  f  E(t)  E(t  +  r)dt 
T-*co  "T 

where  r  is  a  delay  time.  Figure  26  shows  the  autocorrelation  function 
of  the  hot-wire  bridge  voltage  output  for  a  few  centerline  locations. 
Since  <>(7)  is  a  pure  function  of  delay  time,  the  integral  time  scale  of 
the  turbulence  can  be  defined  as  the  area  under  the  autocorrelation 
function  curve  (42).  As  shown  in  the  tabulation  in  this  figure,  the 
integrated  time  is  on  the  order  of  0.1  second.  In  other  words,  the 
characteristic  frequency  of  the  turbulence  in  the  flow  field  of  the 
present  experiment  is  low  (on  the  order  of  10  Hz)  which  is  well  within 
the  capability  of  the  electronic  instruments  used. 


We  conclude  chat  the  turbulence  in  the  primary  flow  region  is 
characterised  by  low  frequency  (~10  Hz)  and  large  amplitude  (~ 1  ft/sec) 
fluctuations. 

Hot-Wire  Response  Equation 

Since  the  cooling  effect  on  a  hot-wire  depends  primarily  on  the 
velocity  component  Wj_  perpendicular  to  the  wire  and  to  much  less  extent 
on  the  parallel  component  ,  it  is  useful  to  define  an  effective 
instantaneous  velocity  W  (44) : 

W2  =  Wj_2  +  k2W„2  (5-1) 

where  k  is  an  experimental  constant.  For  a  wire  of  length- to-diameter 
ratio  greater  than  600,  k2  is  essentially  zero  (44).  Equation  (5-1)  will 
be  used  to  obtain  the  value  of  W  in  Equation  (3-2).  Rearranging 
Equation  (3-2),  the  following  relations  are  obtained: 


2 

=  W  (5-2) 

I 

In  these  equations  the  ambient  tempterature  dependent  quantities  A  and 
B  are  lumped  on  the  left  hand  side.  We  will  assume  and  later  verify  on 
the  basis  of  computer  simulated  output  calculations  that  the  following 
ratios  are  substantially  independent  of  the  fluctuations  in  A  and  B. 
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where  an  over  bar  refers  to  the  corresponding  time-averaged  value. 

are  primarily  functions  of  the  normalized  turbulence  velocity  and 
the  wire  orientation  angle  $  . 

Data  in  Primary  Flow  Region 

Eight  locations  in  the  primary  flow  region  were  probed.  The  data 
consisted  of  the  time-averaged  values  (averaging  time  ~2-5  minutes)  of 
E,  E2,  W  and  W2  for  $w  from  -90°  to  +90°  at  15°  intervals.  To  facilitate 
the  data  taking,  the  temperature  dependent  quantity  A  in  Equations  (5-2) 
was  replaced  by  a  fixed  value  A'  nearly  equal  to  A.  The  data  thus 
obtained  were  used  to  compute  the  correct  values  of  W  and  W2,  The 
conversion  formulas  are  given  below  (see  Appendix  1/  for  details): 


w  ■  ^  w  -  ic,2y112  -  c22 

W2  *  C12W*2  -  6C2W(C2  +  2W1/2)  +  4C2W1/2(2w1/2 2 


(5-4) 
C22)  -  c24 


where  =  B’2,  C2  =  (A-A’)/B,  W  =  (E2-A')2,  W'2  »  (E2-A')4  . 

For  all  the  data  taken,  A'  was  set  equal  to  25.  At  this  A1  setting, 
the  whole  data  taking  system  (Figure  15)  was  calibrated  against  a  known 
d.c.  plus  a.c.  voltage  input.  The  calibrations  of  E2,  W' ,  W'2  were 
plotted  versus  their  corresponding  line  outputs  on  the  recorder  and  are 
shown  in  Figure  27. 
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Since  A'  was  taken  as  a  constant,  we  are  essentially  treating  A 
and  B  as  constants  evaluated  at  the  average  ambient  temperature.  The 
maximum  error  introduced  due  to  this  assumption  may  be  estimated  from 
the  dependences  of  A  and  B  on  fluid  properties  (p.  79  ,  reference  42)  and 
can  be  shown  to  be  on  the  order  of  1%.  As  a  check  the  effects  of  the 
fluctuations  of  A  and  B  due  to  temperature  fluctuations  is  evaluated  by 
numerical  simulations  as  will  be  described  later. 

The  hot-wire  data  obtained  are  shown  in  Figure  28.  The  ratios 
defined  in  Equations  (5-3)  at  *w  =  90°  may  be  used  to  describe  the 
shape  of  these  curves  in  Figure  28.  These  "shape  parameters"  and  the 
relevant  hot-wire  data  are  summarized  in  Table  4.  A  comparison  between 
W  and  <?w  shows  clearly  that  the  turbulence  intensity  is  indeed  much 
larger  than  that  reported  for  typical  low  turbulence  measurements  (e.g., 
references  42,  45). 

Indraft  Data 

Two  indraft  profiles  at  the  edge  of  the  plate  were  measured.  These 
measurements  were  all  made  in  the  secondary  flow  region,  a  laminar  flow 
region.  Therefore,  direct  measurement  using  only  the  hot-wire  bridge 
voltage  output  E  was  possible.  Since  the  velocity  in  this  region  was 
small,  a  relatively  high  wire  operating  temperature  was  desirable.  The 
wire  was  operated  at  approximately  400°F.  A  set  of  calibration  curves 
was  obtained  for  this  purpose  (Figure  13).  Figure  13  shows  that  the 
velocity  vector  at  any  location  can  be  determined  once  the  hot-wire  bridge 
voltage  outputs  for  *w  =  0°  and  90°  (Eh  and  Ey  respectively  are  obtained. 
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The  results  are  shown  in  Figure  29.  Clearly,  the  largest  indraft 
occurs  near  the  ground  level. 

NUMERICAL  SIMULATION  OF  HOT-WIRE  DATA 

k—w— rwia  i  mmmwmm — a— — i a— — — <— i— ttmmmmmm we— I— » 

Described  here  is  an  Idealised  numerical  model  to  generate  simulated 
hot-wire  data  for  a  given  set  of  flow  parameters.  The  results  thus 
obtained  will  be  used  to  interpret  data  from  experimental  measurements. 
The  Numerical  Model 

Consider  the  mean  velocity  vector  Vm,  the  fluctuating  velocities 
(u,  v,  w)  and  the  hot-wire  described  in  the  following  coordinate  system: 


where  v  is  the  random  fluctuating  velocity  component  parallel  to  Vm,  w 
is  the  random  fluctuating  velocity  component  in  the  horizontal  direction 
normal  to  v,  and  u  is  the  random  fluctuating  velocity  component  in  the 
direction  normal  to  both  v  and  w  such  that  (u,  v,  w)  form  a  right-handed 
triad.  The  direction  cosines  of  Vm,  u,  v,  w  and  the  wire  in  the  (x,  y, 
z)  coordinate  system  are  given  below: 


Denoting  the  component  of  u,  Vm+v,  and  w  in  the  direction  of  the  wire, 
respectively,  by  E^,  E2,  and  E3  the  above  tabulation  gives,  after  a 
little  manipulation: 

Ex  *  si~(tw-^)  +  sin*v  cos*w  [i-sinOw+«v)] 

E2  *  c08*w  C08<*w+«v>  (5-5) 

E3  -  cos(*w-^)  -  cos*y  cos#w  [l-sin(*wi-fy)J 


Finally,  for  given  values  of  Vffl*  (u,  v,  w) : 

Dj|  -  [aEj  +  (,  +  Vb>E2  +  wE3]  2 
-  u2  +  (v  +  Vm)2  +  w2  -  wj 


(5-6) 


These  expressions  will  be  used  in  Equations  (5-1)  and  (5-2)  to  gener&ce 
the  desired  hot-wire  output. 

For  convenience  of  specification  of  the  fluctuations,  the  following 
normalized  variables  are  defined: 


9 


(5-8) 
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are  the  variance  of  the  fluctuations.  The  correlation  coefficients  of 
turbulence  may  be  written  in  terms  of  the  above  quantities; 


K. 


uv 


uv  XX  u*f 

— —  *  uv  ,  K  .  *  — — 

9  9  Uff  9  9. 

UV  U  t 


•  u#  ,  etc. 


(5-9) 


The  values  of  these  correlation  coefficients  may  be  positive  or  nega¬ 
tive  but  srnist  not  exceed  unity  in  magnitude. 

Substituting  Equation.!  (5-1) »  (5-6)  and  (5-7)  into  (3-2)  and  making 
modifications  for  studying  the  effects  of  temperature  dependent  fluctua¬ 
tions  of  A  and  B,  we  have: 


E2  »  (1  +  *9f 7)  A  +  (1  +  799^7)  BVm 


1/2 


Wj. 

Vm 


!+  k2  W ' 

\vm/ 


2n  1/4 


(5-10) 


For  computations,  we  can  without  loss  of  generality  take  «  *=  -1  and 

V  =1.  The  factor  7  accounts  for  the  fact  that  A  and  B  do  not  have 
m 

the  same  dependence  on  temperature. 

The  results,  with  combinations  of  the  »'s  and  K’s  as  defined  in 
Equations  (5-8)  and  (5-9),  will  be  used  for  this  study.  The  flow 
variables  characterizing  the  turbulence  (u,  v,  w,  ®')  define  actually 
random  variables,  sequences  of  which  can  be  simulated  with  the  numbers 
ODtained  from  a  random  number  generator  (see  details  in  Appendix  E), 
given  the  values  of  a's  and  K's.  After  making  some  idealizations 
described  below,  the  desired  hot-wire  response  and  its  moments  can  be 
readily  computed  as  functions  of  the  flow  variables  from  the  above 
equations. 
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Soma  Numerical  Results 

Calculations  were  made  for  four  types  of  probability  density  distri¬ 
butions.  S qae  typical  results  are  presented  in  Tables  5  end  6.  It  is 
seen  that  fairly  good  results  in  both  reproducing  input  data  and  errors 
can  be  obtained  from  calculations  using  400  samples.  The  effects  of 
thd  type  of  distribution  function  appear  minor. 

Figures  30  and  31  show  some  simulated  hot-wire  data  corresponding 
to  the  experimental  results  given  in  Figure  28.  Resemblance  in  basic 

shape  of  the  data  curves  in  Figures  28  and  30  and  31  is  clearly  shown. 

2 

For  the  hot-wire  used  in  this  experiment,  the  value  of  k  was 
found  to  be  essentially  equal  to  zero.  For  data  reduction,  the  shape 
parameters  and  the  effective  flow  velocity  for  the  case:  k*  *  0,  zero 
correlation  and  temperature  fluctuations,  Gaussian  distribution,  are 
presented  in  Figure  32. 

DATA  REDUCTION 

Data  reduction  utilizes  the  results  of  simulated  hot-wire  data 
based  on  the  following  idealizations:  k^  ■  0,  no  correlations,  no 
temperature  fluctuations,  and  Gaussian  distribution  of  fluctuation  of 
each  velocity  component.  The  justification  of  these  assumptions  will 
be  discussed  under  the  next  heading. 

The  data  reduction  procedure  is  described  below  with  an  actual 


example. 


Table  5.  Comparison  of  simulated  hot-wire  data  -  basic  output. 


1C  refers  to  the  type  of  random  number  probability  distribution.  1C  *  -1  (normal  distribution  scheme) 
1C  =  0  (general  scheme),  IC  =  1  {Eq.  E-8),  IC  =  2  (Eq.  E-9) ,  IC  “  3  (Eq.  E-10). 


Table  6.  Comparisons  of  simulated  hot-wire  data  -  correlation 
coefficients , 


iO 


lA 

o 

u 

3 

i 

> 

o° 

•  « 

H 

«  •=* 

O  (i 

• 

& 

It  •» 

fi2  * 

>o„> 

II  II 

m 

u 

M  lOO 

J*  * 

O  It 

o  \3 
•  « 

U 

: 

*  - 
*•  *  »  4 

«  ©  O  p  « 
w  •  •  5  «d 

Id  N  H  *Q 

*, a  *  4  5 

g***!  a- 

•  M 

M  M 

CM 
/-s 
<Si  U} 

tz  0- 

Vr^ 

2.927 

3.358 

3.124 

2.440 

3.225 

2.836 

2.612 

2.751 

2.595 

1.442 

5.282 

5.117 

w 

f  ps 

— 

1.296 

1.392 

1.347 

1.313 

1.417 

1.354 

1.320 

1,372 

1.319 

Ov  CA  CO 

OO  O  <T 

©  <a  <a 

•  •  * 

H  H  H 

-o 

o 

t4 

to 

•0 

A 

CM 

\ 

H 

|<M  /-s 

~h*  CO 

*4  a 

1  ^ 
w 

1.053 

1.100 

1.090 

1.101 

1.121 

1.101 

1.083 

1.112 

1.085 

1.019 

1.055 

1.074 

_ 

E 

volt 

1.747 

1.757 

1.754 

1.759 

1.763 

1.758 

1.753 

1.761 

1.753 

1.736 

1.743 

1.749 

I  correlation  coefficients  J 

'"\l 

-.121 

-.109 

-.099 

.060 

-.020 

-.005 

-.051 

-.039 

-.014 

-.010 

.011 

.018 

3 

cm  ia  ^ 

§  2  g 

•  •  • 

SB 

SB 

3 

o  ia  <a 

O  CM  CM 

*0  4A  tA 

•  •  • 

.751 

.496 

.502 

f-4  CNI  VO 

vo  f*  ph 

*  »A  < 

•  •  • 

CO  ^  VO 

*■4  fM  OV 

rn  pH  vo 

• 

j 

Ike 

vMhI 

H| 

kb 

■ 

FffM 

mflA 

j 

.031 

-.015 

-.016 

-.113 

-.014 

.011 

-.013 

-.054 

-.100 

rH  rn  m 

•O  tA  CM 

©  ©  © 

«  •  • 

i 

> 

<f  vO  00 

00  00  04 

lA  iA  lA 

•  •  • 

.481 

.500 

.529 

CM  O  t4 

so  m  »■* 

*A  ia 

•  •  • 

.290 

.256 

.351 

(0 

g 

4J 

s 

4J 

o 

p 

f4 

<U 

to 

8 

b 

IB 

.125 

.128 

.122 

.120 

.123 

.121 

.120 

.120 

.120 

.091 

.145 

.126 

<0 

S  g  o 

lA  lA  kA 

CM  Ph  CM 

00  00  O 

<r  tn 
•  •  • 

•O  ©  o 

O'  O  O 

*0  *A  vA 

•O  vO  0* 

O  CM  pH 

-o 

(0 

a  a. 

%  V* 

.527 

.522 

.506 

.505 

.503 

.506 

oo  oo  ov 

O'  O'  O' 

*0 

•  •  • 

O'  CM  O' 

O'  <r  <f 

CA 

♦  •  * 

IB 

>  O. 

1.054 

1.070 

1.011 

0.990 

1.016 

1.009 

1.000 

0.995 

0.999 

0.735 

1.192 

1.041 

*sidm»s 

JO 

j»qun\ 

100 

400 

1000 

100 

400 

1000 

100 

400 

1000 

100 

400 

1000 

% 

O 

H 

CM 

■ 

See  footnote  of  Table 


51 


1.  Center  line  Locations 

On  the  centerline,  symmetry  requires  that  *  ov  and  U  *  0,  so 
the  shape  parameters  plotted  in  Figure  32  are  directly  applicable. 
Consider  the  data  for  X  *  0,  Y  «  3": 


Shape  Parameter 
of  Actual  Data 

Range  of  <ru/Vm 
from  Figure  32 

Common  Range 
of  »u/Vm 

0.775 

0.52-0.55 

0,52-0.55 

U 

0.52-0.55 

0.608 

0.51-0.56 

The  value  of  <ru  will  be  chosen  (within  the  above  common  range)  such 
that  the  values  of  (*u  -  obtained  from  («ru/Vm,  *0>  and 

(cu/Vm,'®.^)  are  the  same,  this  being  a  requirement.  In  this  manner, 
the  value  of  vu/Vm  and  (wu  -  <*v)/Vm  are  found  equal  to  0.52  and  0, 
respectively. 

From  Figure  32  for  »U/Vm  ?**0.52  and  <*u i -t*v)/Vn  *  0. 

W  =  1.18  Vm 

The  hot-wire  data  in  Figure  28  shows  that: 

W  -  1.58  fps 

Finally,  Vm  =  1.35  fps 

ffu  =  ffv  =  %  =  0.7  fps 

2.  Off-Centerline  Locations. 

For  off-centerline  locations,  Vm  is  making  an  angle  4V  with  the 
centerline.  This  angle  is  s^t  approximately  equal  to  the  phase 
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shift  of  the  peak  value  of  the  hot-wire  data.  Using  the  peak  values 
as  a  basis, &g,  (Sip  and  W/Vffi  are  then  obtained.  The  rest  of  the 
data  reduction  is  carried  out  in  the  same  manner  as  described  in  (1). 

A  summary  of  the  reduced  data  is  given  in  Table  7.  It  is  seen 
that  the  turbulence  velocities  are  indeed  large  (of . comparable  magnitude 
as  the  mean  flow  velocity),  and  that  the  turbulence  decreases,  as  it 
must,  as  the  distance  from  the  plate  increases.  It  is  remarked  that  in 
the  data  reduction  procedure  described  above  the  effects  of  local  free 
convection  due  to  the  hot-wire  has  been  ignored,  i.e.,  the  hot-wire 
output  is  assumed  unaffected  by  the  wire  orientation  with  respect  to 
the  mean  flow.  In  view  of  the  large  turbulence  velocities,  this  appears 
to  be  justified  (see  Chapter  III). 


Table  7,  Sunaaary  of  hot-wire  velocity  results  (  m  520°F). 
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If  we  assume  that  the  turbulence  is  characterized  by  ?u  and  the  time 
obtained  from  the  autocorrelation  function,  then  au^r can  be  used  to 
estimate  the  order  of  magnitude  of  the  eddy  diffusivity.  Table  7  f.hovs 
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chat  the  values  thus  calculated  are  consistent  with  the  estimates  made 
based  primarily  on  temperature  data  (Table  3). 

JUSTIFICATIONS 

Simulated  hot-wire  data  have  been  obtained  for  a  variety  of  cases 
based  on  different  assumptions  such  as :  type  of  random  number  distribu¬ 
tion  function,  correlation  coefficients,  temperature  fluctuations,  and 
the  values  of  A  and  B.  The  results  show  that  the  computed  values  of  «q 
and  12^  are  somewhat  insensitive  to  the  assumptions  used  and  the  values 
of  $2  and  W/Vm  show  more  pronounced  dependence  on  assumptions. 

Figure  33  shows  comparisons  of  some  of  the  typical  results  with 
the  largest  differences  (l.e.,^  and  W/Vm).  Clearly,  relatively  large 
differences  (although  small  by  themselves)  may  be  caused  by  the  type  of 
random  number  distribution  functions  used  and  the  results  are  relatively 
insensitive  to  other  assumptions.  In  the  data  reduction  procedure,  only 
the  results  of  ^  and  U/Vm  are  used.  In  view  of  their  insensitivity 
to  the  assumptions,  this  procedure  is  reasonably  general  and  the  expected 
error  involved  due  to  assumptions  may  be  best  shown  by  that  of  Vm. 

In  order  to  show  quantitatively  the  effects  of  assumptions  on  the 
reduced  data,  the  procedure  for  the  actual  data  reduction  was  repeated 
for  a  number  of  cases  using  simulated  hot-wire  data  based  on  different 
assumptions.  Relative  errors  are  defined  as  the  percentage  deviation  of 
the  ct' -’listed  results  from  their  corresponding  values  tabulated  in 
Table  7.  Results  show  that  the  relative  errors  of  the  reduced  turbulence 


54 


velocities  0^,  *v  and  the  mean  flow  velocity  Vn  are  of  comparable 

magnitudes.  This  is  expected  because  of  the  dominant  dependence  of 

errors  on  W/\L  described  above, 
m 

Table  8  shows  the  relative  errors  of  V  and  one  case  of  0,.  It  can 

m  u 

be  seen  that  the  possible  error  in  data  reduction  due  to  the  various 
assumptions  used  is  on  the  order  of  10%.  The  mean  flow  direction  is 
unaffected  by  the  assumptions  in  any  way  because  it  is  determined  by 
the  location  of  the  peak  value  of  the  actual  hot-wire  data.  The  error 
of  the  mean  flow  direction  is  mechanical  which  is  less  than  ±2°. 


VI.  NUMERICAL  MODEL 


The  purpose  o£  this  chapter  is  to  present  e  numerical  method  for 
predicting  the  flow  field  near  the  plate. 

MATHEMATICAL  FORMULATION 

In  the  cylindrical  coordinate  system  sketched  below, 
y.v 


the  flow  field  to  be  described  is  that  generated  in  the  half  space 
y  ;»  0  by  a  heat  source  of  diameter  D  at  y  •  0.  An  order  of  magnitude 
analysis  shows  that  the  viscous  dissipation  (Eckert  number  M.0"®  for 
our  problem)  and  the  flow  work  terms  in  the  energy  equation  are  negligible 
(Chapter  14,  Reference  45).  In  the  primary  flow  region,  the  region  in 
which  there  is  a  dominant  upward  moving  stream  of  eddies,  the  density 
(or  absolute  temperature)  variations  are  small  compared  with  the  ambient 
density,  therefore,  the  familiar  Boussinesq  approximations  are  used  (1). 

At  present,  the  physical  behavior  of  the  flow  field  is  not  well-under- 
stood,  nor  have  there  been  any  proposed  turbulent  transport  formulations 
in  the  literature.  The  constant  eddy  diffusivity  model  is  the  simplest 
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plausible  turbulent  transport  model  and  is  chosen  for  the  calculations 
described  below.  After  introducing  these  simplifications,  and  a  charac¬ 
teristic  velocity  Uc^  *  the  resulting  equations  for  the  mean  flow 

field  become  the  same  as  those  for  laminar  free  convection  (13,  36). 
These  differential  equations  in  dimensionless  form  are  given  below: 
Continuity 

+  (6-D 


Momentum 


8u  t  du  du  d  p  ,  1  /d^u  1  du  u  d^u 

dv  ,  dv  dv  dp,  1  cfiv  ,  1  dv  ,  d?v\  ,  _ 

^  +  u^+vd7=“^  +  vg?;  fc  +  x^  +  a5?J+T 


(6-2) 

(6-3) 


Energy 


dT  dT  ,  dT  1  Id2 T  1  dT  A  d2T 

5T  +  u  K  +  v  37  *  pTv^r  ssr +  x  5T  + 


(6-4) 


where  x,  y,  u,  v,  p,  T,  GrT,  Pr  and  t  are,  respectively,  the  dimensionless 

radial  and  axial  coordinate  (x  =  X/D,  y  =  Y/D),  dimensionless  x  and  y 

component  of  the  mean  flow  velocity  (u  =  U/U„,  v  =  V/U„),  pressure 

c  c 

2 

coefficient  (p  *  (P-P^/Poo^c  )>  dimensionless  temperature  (T  «=  (P-^/AP-j.), 
turbulent  Grashof  number  (GrT  *  g/SAft^D-fy^),  turbulent  Prandtl  number 
(Pr  =  »  a°d  dimensionless  time  (t  =  rUc/D).  AAp,  and  are, 

respectively,  the  equivalent  temperature  rise,  momentum  and  thermal 
dirfusivities  for  turbulent  flow. 

The  physical  boundary  conditions  are:  Zero  disturbance  at  large 
distance  from  the  heat  source  and  non-slip  and  thermally  insulated 
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conditions  at  the  ground  surface  (y  *  0),  and  T  •  i  at  the  source. 


METHOD  OF  SOLUTION 

Defining  a  stream  function  S  such  that: 

as 


xu 


-XV 


ay’ 

as 

3* 


(6-5) 


automatically  satisfies  Equation  (6-1).  Subtracting  the  derivative  of 
Equation  (6-3)  with  respect  to  x  from  the  derivative  of  Equation  (6-2) 
with  respect  to  y  and  making  use  of  Equation  (6-1)  yields: 


az  .  az  .  az 
3_+us.+v- 


3y  x 

The  new  quantity  Z,  defined  as: 


i  ax  . 


r*2, 


z  3  az 

+  x  Sx  + 


_  du  dv 

“xZ  •  SF  ‘  ^ 


(6-6) 


(6-7) 


in  Equation  (6-6)  is  called  the  modified  vorticity  function. 
Substituting  Equation  (6-5)  into  Equation  (6-7)  a  relation  between 
the  stream  function  and  the  modified  vorticity  function  is  obtained: 

-x2Z  « 


is  i  as  ,  a2s 

-  X  5T  + 


(6-8) 


The  set  of  auxiliary  conditions  for  the  present  problem  becomes: 

(a) t<0:  u«v-T*S»Z*0  for  y  0. 

(b)  t  £  0: 

u,  v,  T,  S,  Z  -  0  at  large  distances  from  the  source, 
y  «  0: 

T  -  1  |x|  <  D/2 

^  “  0  »xi  >  D/2 


>  (6-9) 
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u  =  v  =  S  =  0 


u 

Z 


S  "fe  ° 


(6-9) 


) 


The  set  of  equations  to  be  solved  consists  of  (6-4),  (6-6),  (6-8)  and 
(6-5).  These  equations  were  solved  numerically. 

We  seek  steady  state  solutions  for  the  problem  posed  above.  In 
general,  two  basically  different  numerical  approaches  may  be  considered 
the  time- independent  and  the  time-dependent  formulations  and  solutions. 
The  former  approach  requires  the  simultaneous  solution  of  the  steady 
state  version  of  all  the  governing  equations.  This  requires  iterative 
procedure  which  may  be  complex  and  time-consuming  (13,  35,  36).  The 
latter  approach  is  to  introduce  a  disturbance  (i.e.,  for  our  problem, 
temperature  at  the  source)  into  the  initially  undisturbed  region  of 
interest  at  t  >  0.  The  subsequent  changes  in  the  region  of  interest 
as  the  time  increases  will  be  described  by  the  governing  equations. 

The  steady  state  solution,  if  any,  will  be  that  obtained  at  large  time. 

Various  types  of  number ical  flow  calculations  (for  example, 
references  9,  34,  43,  46,  47,  48,  49,  etc.)  have  shown  that  using  a 
time-dependent  numerical  scheme  to  obtain  a  steady  state  solution  is 
convenient.  The  time -dependent  approach  was  used  here. 
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FINITE  DIFFERENCE  SCHEME 
Grid  System 

For  an  axially  symmetric  problem,  calculations  need  be  performed 
only  over  a  radial  plane  containing  the  symmetry  axis.  The  square  mesh 
grid  system  for  such  numerical  computation  is  described  below: 


y,v 

4 


j  *  jmax 


j  -  1 


axis  of  symmetry 
'(or  centerline) 

^ upper  boundary 

-  -  --n 

dx  ”  Ay 


-+\  Ay  U- 

i  f  *”t^4-(i,j+i) 

Tax  I  j  J 

l-l-l 


source 


; 


ground  surface 


% 


outer 

boundary 


77 T7 T7TT7 7^'V  7 x,u 

(radial 

i  ■  ®  i  “  imax  axis) 


Henceforth,  two  subscripts  (i,j)  will  be  attached  to  each  dependent 
variable  (when  necessary)  for  identification  purposes.  Thus  T^j  refers 
to  a  dimensionless  temperature  on  the  line  x  -  0. 

For  identifying  the  particular  size  of  array  of  nodal  points  and 
the  relative  size  of  the  source,  the  following  symbol  will  be  used  in 


the  future  discussions: 


V**"n5P 


61 


"ioax  x  jmax,  1-m" 

Here  the  first  part  specifies  the  array  size  and  the  second  part 
means  that  the  source  occupies  the  region  from  i  *  1  to  i  *  m,  or  source 
size  *  (m-l)dx. 

Difference  Forms 

The  derivatives  in  the  differential  equations  to  be  solved  will  all 
be  approximated  by  their  difference  analogs.  Three  difference  forms 
(forward ,  backward,  and  central)  will  be  usedofov  the  first  derivatives 
and  a  single  form  will  be  used  for  the  second  derivatives.  To  illustrate, 
the  space  derivative  for  any  function  S(x,  y,  t)  are  given  below: 

(a)  First  derivatives 


Forward  difference 


Backward  difference 
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Stability  Criterion  and  the  Difference  Equations 

Equations  (6-4),  (6-6),  (6-8),  and  (6-5)  are  replaced  by  their 
difference  analogs,  Wilkes  (9)  has  investigated  this  system  of  equations 
extensively  for  a  free  convection  problem  in  a  rectangular  cavity  with 
heated  vertical  walls.  He  used  a  differencing  scheme  (50)  such  that, 
for  the  convective  terms  on  the  left  hand  side  of  Equations  (6-4)  and 
(6-6),  the  forward  difference  form  was  used  when  the  coefficient  velocity 
was  negative,  and  the  backward  difference  form  was  used  when  the 
coefficient  velocity  was  positive.  The  central  difference  form  was  used 
for  all  other  terms.  Based  on  a  linear  Fourier  analysis  the  following 
stability  criterion  may  be  derived: 

ttS[%i+lwL+_^.(a7+a7]j'1  (6-io) 

Taking  u  as  negative  and  v  as  positive,  for  illustrative  purposes,  the 
difference  forms  of  Equations  (6-4)  and  (6-6)  are  written  as  follows: 


63 


dZ 

3t 


-  Z4 


Zj  i  -  Zj 


=  -  u 


i.j 

+ 


i,j  Ax 


2xAx 


2xAx 


1 

'  i.j  Ay 

Zi+l,j  “  2Zi,j  +  zi-l,j 

v£'rT 

Ax2 

~  *i 

i-lJ 

,  zi,j-n •  2zi,i +  zi,j-i. 

T  lyl 

(6-12) 


Equation  (6-8)  may  be  written  first  as  follows: 


S...  i  ~  2S.  .  +  S.  .  ,  S.,.  .  -  S.  ,  . 

^l  >  j  l»j  l“l>j  1*^1 1 J  r-l,j 

- ^2 - 2xAx 

S,  .  .  -  2S.  .  +  S.  ,  .  - 

+  *  -  X2  Z. 

Ay2  i,j 

In  the  above  equation,  ^  wae  solved  by  an  iterative  procedure  using 

an  over-relaxation  scheme  (51).  Rearranging  and  introducing  the  over¬ 
relaxation  parameter  W,  and  x  *  (i-l)Ax,  we  have: 


Here,  the  left  hand  side  is  the  value  of  ^  from  the  most  recent 

iteration  based  on  the  established  values  of  S?  etc.  from  the 

1 » J 

previous  iteration  on  the  right  hand  side.  Finally,  from  equation 
(6-5): 
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u 


.  Si.W  •  sl,i-l 
i» j  2x4y 


(6-14) 


i,j  2x4x 


Boundary  Conditions 

The  governing  differential  equations  for  this  problem  are  of  the 
elliptic  type;  the  specification  of  conditions  on  all  the  boundaries 
is  required.  We  define  a  natural  boundary  as  one  at  which  appropriate 
boundary  conditions  are  obvious.  The  only  natural  and  finite  boundaries 
in  this  problem  are  the  ground  surface  and  the  axis  of  symmetry.  The 
others  are  at  infinity.  For  numerical  calculations,  one  is  restricted 
to  perform  calculations  in  a  finite  domain  of  a  feasible  number  of  nodal 
points.  To  meet  this  restriction  one  may  either  map  the  semi-infinite 
space  y  >  0,  into  a  region  defined  between  0  and  1  or  assume  some 
approximate  artificial  boundaries.  Mathematically  it  is  possible  to 
map  the  infinite  domain  into  a  finite  one  by  a  transformation  such  as 
the  following  pair: 


i 


x 


1  +  C.  x  * 


ci  i  -  r 


0  £  X  0£  f  Si. 


c2y 


(6-15) 


1  +  C2  y  » 


c2  1  -  1 


;  o  £  y  £  00,  0  2s  *  £l. 


where  and  C2  are  constants  which  determine  respectively  the  manner 
in  which  the  x-  and  y-coordinate  are  to  be  shrunk. 
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Some  attempts  were  made  to  use  the  transformed  coordinate  system  as 
described  above.  It  was  found  that  although  the  calculation  seemed  to 
converge  reasonably  well  and  produce  stationary  results  near  the  source, 
it  exhibited  large  oscillations  at  large  distances  from  the  source. 

This  is  because  of  the  singular  behavior  near  t  ■  1  and  v  -  1  and  the 
associated  unusually  large  truncation  error  due  to  the  transformation 
in  these  regions.  Figure  34  is  a  plot  of  some  numerical  results  using 
the  mapped  coordinates  for  the  case  of  Gr-p  *  7.4  x  10^  for  the  different 
stages  of  calculations.  Although  the  Grashof  number  used  in  this 
example  is  too  large  to  be  realistic  in  the  physical  sense,  the  plotted 
results  do  show  the  large  oscillations  near  infinity.  Furthermore,  it 
took  about  50%  more  computing  time  than  when  an  ordinary  physical 
coordinate  system  was  used.  For  these  reasons  the  physical  coordinate 
system  was  used  and  non-natural  boundaries  defined  where  necessary. 

Some  reasonable  boundary  conditions  on  these  non-natural  boundaries 
are  constructed.  Boundary  conditions,  with  discussions  where  appropriate 
are  given  below: 

(a)  Source  Plane  (j=l,  y=0) 

T,  =  1  for  x  <  D/2 

i,l 

=  0  for  x  >  D/2.  Due  to  the  inflow  of  air  at  ambient 

M 

temperature,  this  condition  is  nearly  equivalent  to  T  *  0. 
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Since  no  flow  is  crossing  the  source  plane,  y  *  0  is  a  stream  surface. 
Thus  we  can  set: 


By  Equation  (6-8): 


Expand  S  near  y  *  0, 
Si,2"  Si,l+  *y 


Z.  , 


S.  ,  =  0  . 


1  d2S 


X7  dy2  j 


i»l 


as 

dy 


i.l 


+  iid 
+  2 


a2s 


si,3  “  st,i  + 


Ids 


+ 


4av2  /a2s 


.  ^2  /a3sl 

i,l+  *  M 

i£sl  .  8a£  f£s| 

ffil  i  L  +  6  [dy7/ 1  1 


Noting  that: 


si,i  ’ 0  md  <xu)i,i  *  °- 


ld3s\ 

and  eliminating  [dyTj  fro®  the  above  two  expansions,  we  get: 

Wl 


'i»l 


8si,2  "  si,3 

-  'iiy’i" 


Finally, 


8s‘.,2  '  Sj.3 


2x^Ay^ 

This  expression  (9)  is  accurate  to  the  second  order  of  Ay.  This  high 
accuracy  is  desirable  in  order  to  minimize  the  inaccuracy  in  u  and  v 
as  calculated  through  Equation  (6-5). 

(b)  Axis  of  Symmetry  (i*l,  x=0) 

Based  on  the  symmetry  of  the  flow  field  and  the  requirement 
that  field  variables  and  their  gradients  be  everywhere  continuous, 
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Timax,j  "  0  ‘ 

This  condition  simply  means  that  the  temperature  reaches  the  ambient 
value  on  the  outer  boundary. 

We  have  so  far  specified  zero  temperature  gradients  on  the  other 
three  boundaries  (except  at  the  source).  In  a  steady  flow  problem,  the 
velocity  field  is  established  and  the  energy  equation  becomes  linear. 
The  specification  of  zero  temperature  gradient  on  this  boundary  would 
result  in  a  trivial  solution  of  T  *  1  everywhere  inside  the  boundaries. 
Furthermore,  the  temperature  on  the  outer  boundary  can  be  higher  than 
that  of  the  ambient  only  by  heat  diffusion  in  an  upstream  direction. 


The  characteristic  diffusion  length  is  proportional  to 


U hD\1/2 
\”c/  ’ 


our  problem,  this  diffusion  length  is  on  the  order  of  10"  1  foot. 
Therefore,  at  a  distance  not  far  from  the  plate,  the  air  temperature 
would  be  very  near  the  ambient  value.  Thus  the  choice  of  Tj^]MX  j  “  0 
at  sufficiently  large  value  of  imax  (e.g.,  1.5D)  is  well  justified. 

For  velocity  conditions  at  this  boundary,  we  choose: 

^imax,j  ^imax-l,j 

Wj-°- 

The  first  of  these  implies  v  «  0  and  xu  *  constant  from  Equations  (6>5) 
and  (6-1).  Since  the  largest  velocity  change  takes  place  within  the 
boundary  layer  near  the  ground  surface,  the  remaining  portion  of  the 
secondary  flow  region  at  this  boundary  may  be  considered  approximately 
irrotational.  Experience  in  numerical  calculations  shows  that  the 
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specification  of  Zioax  j  is  insensitive  to  the  solution. 

Following  the  basic  scheme  discussed  in  (a)  above,  the  boundary 

values  for  the  condition  of  vanishing  first  derivative  were  established 

by  3-point  extrapolation  scheme  (53).  Considering  T.  .  for  example,  a 

•*■>3 

Taylo"  series  expansion  gives: 


T0  .  cs  T,  .  +  Ax 

2,J  l»j 


tel  +  +  ^  ii 


i,j  6  Mi(j 


,  ,»T.  ,  +  2te  tel  +!*£  &]  f£! 

1,J  ox  2  1  Ox*  I  6  dx-3/,  , 

'l,j  'l,j 


\i  “  Ti.j +  34,1 


9 Ax2  Id2 T 


+  —  ter  ,  ,  +  6 


27Ax3 


dr 

Since  ^ 
to  give: 


0,  the  above  equations  may  be  solved  simultaneously 


T  »  (1ST  .  -  9T+  2T,  )/ll. 

* » J  *»J  •♦ » J 

Since  the  profiles  of  interest  have  no  rapid  changes  in  slope  along 
(l,j),  this  method  was  found  to  be  very  satisfactory.  Similarly, 

Hmav  be  obtained  for  specifying  Z,  ,.  Since  v  is  a  quantity 

1,J  * 

derived  from  the  stream  function  S  which  in  turn  is  obtained  by  solving 
Equation  (6-8)  from  the  known  values  of  Z,  it  is  more  accurate  to 
compute  Z.  .  directly  by  extrapolation  from  Z,  . ,  Z,  .,  and  Z,  .  This 
scheme  was  used  and  was  found  to  be  both  convenient  and  satisfactory. 

It  is  realized  that  the  conditions  imposed  on  the  outer  and  upper 
boundaries  are  somewhat  artificial.  Further  justifications  will  be  given 
when  the  domain  of  calculations  is  discussed. 
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Calculation  Procedure 

The  governing  equations  were  programmed  in  FORTRAN  IV  computer 
language  and  the  calculations  were  done  on  an  X2M  7094  digital  computer. 
A  block  diagram  describing  the  calculation  procedure  with  discussions 
where  necessary  is  given  in  Appendix  F.  The  computer  program  listing 
is  given  in  Appendix  G. 

VALIDATING  OF  NUMERICAL  METHOD 
The  Diffusion  Equation  Limit 

We  attempt  to  validate  the  numerical  results  by  going  through  a 
limiting  procedure.  As  a  preliminary  test,  the  velocity  field  in  the 
computer  program  was  suppressed  so  that  the  numerical  results  would 
correspond  to  the  solution  for  a  diffusion  problem  with  the  same 
boundary  conditions.  This  procedure  was  used  also  by  Kane  and  Yang  (36) 
in  a  similar  calculation.  An  exact  solution  for  this  problem  in  closed 
form  is  given  below  (54) : 

2  .  -1  f  2 

x"*8xn  7”?  7-J'  1 T'T  '  '"I  * 

L  V(x-l)  +  y  +  V(**l)  +  y 

Three  numerical  results  are  presented  in  Figure  35.  It  is  seen  that  the 
numerical  results  tend  toward  the  exact  solution  as  the  domain  of 
calculation  is  extended.  The  discrepancy  between  numerical  and  exact 
results  is  apparently  due  to  the  artificial  boundary  condition  on  the 
outer  boundary.  This  is  further  illustrated  for  a  case  in  which  the 
lower  portion  of  the  outer  boundary  is  insulated.  The  results  are  also 
shown  in  Figure  35.  With  respect  to  the  thermal  boundary  condition. 
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the  non-flow  example  is  a  severe  test.  An  inflow  of  air  at  ambient 
temperature  tends  to  make  the  condition  *  0  much  more  realistic. 

As  an  additional  test  of  the  computer  program,  solutions  were 
obtained  for  various  specified  degrees  of  flow  activity.  As  discussed 
earlier,  Gr^  is  the  governing  parameter  in  this  problem.  The  flow 
velocities  are  large  for  large  Gr^.  The  numerical  results  for  two 
Grashof  numbers  (10  and  10J)  are  presented  in  Figure  36.  As  must  be 
the  case,  the  solution  of  the  diffusion  equation  (the  case  of  Gr—*0) 
is  again  approached  as  the  Grashof  number  decreases.  At  large  Gr^,,  the 
region  where  the  temperature  differs  significantly  from  the  ambient 
value  is  squeezed  toward  the  centerline.  This  fact  is  in  agreement 
with  experimental  evidence  (Chapter  IV)  and  theory  (5,38).  It  is  also 
indicative  of  the  unimportance  of  the  temperature  specification  on  the 
outer  boundary  as  the  Grashof  number  increases.  We  are  interested  in 
the  Grashof  number  in  the  range  around  10^. 

Effects  of  the  Domain  of  Calculation 

The  required  computing  time  for  obtaining  a  converged  solution  for 
a  given  problem  depends  on  the  size  of  the  array  of  mesh  points. 
Numerical  experiments  were  carried  out  to  test  the  minimum  array  size 


Kane  and  Yang  (36)  studied  the  laminar  free  convection  problem 
due  to  a  heated  horizontal  disk  in  a  full  space  for  very  small 
Gr  (<100)  and  found  that  their  solutions  tend  to  coincide  with 
the  diffusion  equation  solution  as  Gr  decreases,  as  is  true 
with  the  present  solution.  Direct  comparison  of  these  results 
was  not  made  due  to  difference  in  geometrical  configurations. 
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which  would  produce  reasonable  results  for  the  region  of  interest. 

4 

Calculations  were  made  for  the  case  of  Gr^  *  10  with  a  source  radius 
of  4Ax.  Figures  37  and  38  show,  respectively,  the  effects  of  vertical 
extent  and  horizontal  extent  of  calculation  on  the  numerical  results. 

It  only  the  region  defined  by  x  £  1,  y  5s  1  is  of  concern,  these  figures 
show  no  appreciable  difference  on  the  numerical  results  when  the  source 
radius  is  4ax  and  array  size  is  Increased  beyond  15  x  17.  For  this 
reason  the  ”15  x  17,  1-5"  (see  p.  61  for  explanation  of  this  identification) 
array  was  used  for  further  calculations. 

Effect  of  Mesh  Spacing 

The  truncation  error  for  the  present  system  of  equations  is  on  the 
order  of  K^At  4-  l^Ax  where  and  ^  are  positive  constants  depending 
on  the  dependent  variables  T  and  Z.  By  making  either  At  er  Ax  or  both 
small  the  accuracy  of  the  computed  results  may  be  improved.  Since  we 
are  interested  in  the  steady  state  solution,  the  above  dependence  shows 
that  reducing  Ax  is  an  effective  method  of  improving  accuracy.  A  brief 
Investigation  of  the  effects  of  Ax  on  the  numerical  results  is  described 
below. 

Calculations  were  made  for  GrT  ■  10^  and  10^  by  subdividing  the 
mesh  spacing  in  steps;  the  converged  solution  for  the  "15  x  17,  1-5" 
grid  system  was  obtained  first.  This  first  solution  was  then  used  as 
the  input  for  a  new  system  with  a  mesh  spacing  equal  to  one-half  that 
of  the  first  one,  etc.  In  this  fashion,  the  mesh  spacing  could  be 
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reduced  to  a  reasonably  small  value  with  a  reasonable  amount  of  computing 
time.  A  suumary  of  the  computing  time  for  various  grid  systems  is  given 
in  Appendix  F. 

Figure  39  shows  the  centerline  temperature  distribution  for  Cr^,  *  10 
for  three  values  of  Ax  (3",  1.5",  0.75").  A  nearly  linear  variation 
for  T  with  respect  to  Ax  is  shown  in  this  figure.  The  line  for  Ax  *  0 
was  obtained  by  constructing  a  smooth  curve  for  T(y)  versus  Ax  on  a 
sheet  of  graph  paper  and  reading  the  intersecting  point  of  this  curve 
and  the  line  Ax  ■  0. 

Heat  Transfer 

The  mathematical  formulation  for  our  numerical  calculations  is 
essentially  that  for  a  laminar  flow  case.  Therefore,  it  is  well  to  use 
selected  laminar  flow  results  to  further  check  the  correctness  of  our 
numerical  scheme.  Letting  A •  be  the  difference  between  the  mean  plate 
temperature  and  that  of  the  ambient,  we  have: 


where  h,  k,  6  and  Y  are,  respectively,  the  average  heat  transfer 
coefficient,  thermal  conductivity,  temperature  and  distance  from  the 
source  surface.  If  we  define,  as  before,  T  =  (9  -  6W)  / A  9,  y  =  V/D  and 
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where  T  la  the  dimensionless  mean  temperature  gradient  at  the  source 

y 

surface.  Thus,  the  magnitude  of  the  dimensionless  mean  temperature 
gradient  is  actually  the  Nusselt  number.  Free  convection  heat  transfer 
is  usually  correlated  in  the  following  form: 

Nu  -  C*Grn  or  -Ty  -  C-Gr° 

For  laminar  flow  n  ■  0.25  (e.g.,  reference  27).  This  will  be  used  as 
our  basis  for  comparison. 

Calculations  were  made  for  a  range  of  Gr^  using  the  "15  x  17,  1-5" 

grid  system.  The  variation  of  with  GrT  is  shown  in  Figure  40.  Since 

large  Gr^  also  means  large  heat  transfer,  the  value  of  Ty  as  shown  in 
Figure  40  increases  with  Gr^  as  expected  but  the  value  of  n  was  found 
to  be  0.11.  The  effect  of  mesh  spacing  on  Ty  is  shown  in  Figure  41. 

It  is  clear  that  there  is  an  almost  linear  variation  of  with  respect 
to  Ax.  If  the  values  of  (IL..)^  (see  Figure  41)  are  used,  n  »  0.24, 

C  ■  0.51  may  be  calculated.  This  value  of  n  agrees  well  with  values 

reported  in  the  literature.  No  value  of  C  for  a  horizontal  plate 

(C  ■  0.36  for  vertical  plate)  has  been  reported  in  the  literature,  but 
the  present  result  appears  to  be  good  to  within  a  proportionality 
constant.  For  turbulent  free  convection  n  ■  1/3,  therefore,  some 
discrepancy  between  numerical  results  and  experimental  data  would  be 
expected. 
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NUMERICAL  RESULTS 

Velocity  and  Temperature  Profiles 

Based  on  the  numerical  results  for  three  mesh  spacings  and  the 
graphical  extrapolation  method  described  above,  the  temperature  and 
velocity  fields  for  Gr^,  *  10^  and  10^  extrapolated  to  Ax  *  0  were 
obtained.  These  results  are  shown  in  Figures  42  and  43.  As  expected, 
the  indraft  velocity  is  large  near  the  ground  level  and  the  maximum 
value  of  -u  occurs  always  within  the  plate  radius.  The  outer  edges  of 
the  velocity  and  temperature  profiles  are  more  constricted  toward  the 
centerline  for  large  Gr^  than  small  Gr^,. 

Indraft  Calculations  in  the  Region  near  the  Ground  (l<j<2) 

The  indraft  velocity  near  the  ground  is  ->n *  of  the  most  important 
features  of  the  flow  field.  Measurements  show  that  the  indraft  is 
largest  near  the  ground  and  decreases  rapidly  as  the  distance  from  the 
ground  increases  (Figure  29).  Due  to  this  boundary  layer  behavior 
near  the  ground,  no  feasible  mesh  spacing  for  numerical  calculations 
(e.g.,  D/16)  is  fine  enough  for  predicting  the  detail  in  this  region. 
In  order  to  obtain  a  complete  flow  field  description,  the  present 
method  must  be  supplemented  by  appropriate  auxiliary  calculations. 

If  we  use  the  numerical  results  obtained  with  a  reasonably  coarse 
mesh  spacing  as  an  indication  of  the  outer  flow  field,  this  indraft 
velocity  may  be  calculated  approximately  by  an  integral  method. 


Here  we  refer  to  the  condition  that  only  the  velocity  gradient 
normal  to  the  plate  is  significant. 
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Consider  the  control  volume  described  below: 


By  the  principle  of  conservation  of  energy  and  using  the  perfect 

p 

gas  law  P*  «  -  constant,  the  following  relation  in  dimensionless 

form  may  be  obtained  for  the  control  volume: 


where 


(6-16) 


and  q,  r,  C^,  P  and  fc  are,  respectively,  the  average  heat  flux,  dimen¬ 
sionless  radius  (r  *  R/D),  specific  heat  at  constant  pressure,  pressure 

and  gas  constant.  In  this  equation  Q  may  be  regarded  as  a  known  quantity 

'^The  conductive  heat  transfer  term  is  neglected  due  to  the  small 

local  temperature  gradients.  The  kinetic  and  potential  energy 
terms  are  very  small  and  are  also  neglected. 
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in  which  the  integral  term  may  be  evaluated  from  numerical  results  and 
the  last  term  may  be  obtained  from  experimental  heat  transfer  data. 

The  function  u(y)  is  found  to  be  such  that  the  equality  of  Equation 
(6-16)  holds. 

The  numerical  results  obtained  with  Ax  extrapolated  to  zero  will 

be  used  in  the  calculations.  The  function  u(y)  for  the  interval 

0  £  y  <  ^  wiH  be  determined  based  on  as  much  known  information  as 

possible.  From  numerical  results,  the  end  values  u  ,  and  u  _  are 

m,l  m,2 

fixed.  For  a  better  approximation,  ~  is  also  to  be  satisfied. 


idyjm,2 


This  may  be  computed  from  u  .  in  a  fashion  similar  to  that  given  in 

m,  j 

the  discussion  of  boundary  conditions.  The  resulting  equation  is: 


TT~  (Hu  o  "  18um  *  +  9u„  ,  -  2u 
6Ay  m,  2  m,J  m,4  ul,5' 


In  the  secondary  flow  region,  because  the  velocities  are  small  and  the 
flow  is  in  a  direction  more  or  less  parallel  to  the  plate  surface,  we 
assume  that  the  boundary  layer  equation  for  flows  with  negligible  pressure 

gradient  is  applicable  in  the  neighborhood  of  (m,l).  This  leads  to  the 

12 

r-y  (reference  45).  Considering  a  fourth 
,  *  m,l 

order  polynomial  in  y  with  vanishing  u  and  its  second  derivative  at 
y  *  0,  the  general  form  of  the  polynomial  may  be  written  as: 


u(y)  =  by  +  dy^  +  ey^  . 


(6-17) 


The  coefficients  b,  d,  and  e  may  be  determined  from  the  known  values 


of  and  Q*  An  Ration  similar  to  (6-16)  may  be  obtained 

from  the  principle  of  conservation  of  mass.  This  may  be  used  to  chec' 
the  results  of  Equation  (6-16). 
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The  calculated  results  for  the  case  of  *  S20°F  Is  shown  in 
Figure  44.  A  discussion  will  be  given  at  the  end  of  Chapter  VII. 
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Equations  (7-1)  and  (7-2)  may  then  be  used  to  solve  for  Al^  and  in 


terms  of  Grt,-T- .andiq. . oihewf o|lewing i ralaetene^aaSnobtained : 

1/3 


q  gr  12/3  Kl 

.<>CP^("VJ  L^J 


Zl  n 

g/»  Pr  q  D* 

*M  "  [(pC)*  (-Ty)GrT. 


1/3 


(7-3) 

(7-4) 


In  the  numerical  calculations,  Ty  was  obtained  as  a  function  of  GrT 
(Figure  40).  Equations  (7-3)  and  (7-4)  show  that  if  q  is  given,  A*j 
and  become  functions  only  of  Gr«j. 

COMPARISONS 

Substituting  the  experimental  value  of  q  into  Equation  (7-3),  the 
values  of  A#T,  and  hence  Uc  may  be  calculated  for  a  number  of  Gr^'s. 

The  A&p  and  Uc  thus  obtained  will  be  used  for  converting  the  dimensionless 
numerical  results  into  physical  quantities.  In  this  manner,  the  experi¬ 
mental  and  numerical  centerline  temperatures  were  first  compared  to 
obtain  the  suitable  values  of  Gr^.  Figure  17  shows  that  the  numberical 
results  for  Gr^  *  40000  and  25000  could  fit  the  experimental  data 
reasonably  well  for  «p  ■  520°F  and  410°F,  respectively.  The  values  of 
Ty,  AAp,  Uc  and  corresponding  to  these  Gr^'s  are  tabulated  in  Table  9. 
It  is  seen  that  the  experimental  and  numerical  values  of  are  in  good 


agreement 


VII.  COMPARISON  OF  NUMERICAL  AND  EXPERIMENTAL  RESULTS 


GENERAL  DISCUSSIONS 

In  the  discussion  in  the  early  chapters,  we  have  shown  that,  in 
the  region  sufficiently  close  to  the  plate  surface  (<0.02",  say),  the 
local  air  temperature  can  be  calculated  based  on  molecular  conduction 
(Figures  21  and  22).  Away  from  this  thin  region,  the  flow  is  quite 
turbulent  (Figure  4)  and  the  flow  field,  therefore,  cannot  be  suitably 
described  by  the  molecular  properties  (i.e.,  Gr) .  It  is  not  a  simple 
task  to  attempt  to  describe  the  detailed  flow  behavior  in  the  neighbor¬ 
hood  of  the  source  by  a  single  overall  transport  formulation. 

In  order  to  gain  some  insight  to  the  basic  features  of  the 
numerical  results  using  the  turbulent  Grashof  number  as  the  governing 
parameter,  a  comparison  is  made  between  the  experimental  data  and  the 
numerical  results  presented  in  Chapters  IV,  V  and  VI,  respectively. 

The  turbulent  Grashof  number  is  defined  as : 


g0A0xD“ 


(7-1) 


where  A&j,  and  are,  respectively,  the  temperature  difference  and  eddy 
diffusivity  characterizing  the  flow  field.  The  mean  heat  flux  may  be 


expressed  in  terms  of  a  mean  temperature  gradient  Ty  through  the 
definition  of  a  thermal  diffusivity  as  follows: 

-  *B  M  -  -  ('V«  57  T1 


q 


(7-2) 
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Table  9.  Comparisons  of  experimental  &  numerical  results. 


Experimental  Data 

Numerical  Results  (Ax 

H 

q 

B/hr-ft2 

«M 

ft2/sec 

- - - 

Gr^ 

»T 

y 

°F 

«c 

fps 

ft2 /sec 

520 

see 

.031* 

40000 

6.27 

86 

3.2 

.0324 

410 

560 

25000 

5.64 

61 

2.8 

.0348 

it 

Average  value  of  «^Fr  in  Table  3. 


Having  established  the  values  of  a#x»  the  comparisons  between  the 
horizontal  temperature  profiles  and  updraft  velocities  are  then  compared. 
These  are  shown  in  Figures  18,  19  and  45.  The  numerical  results  show, 
in  general,  more  gradual  profiles  than  do  the  experimental  data.  Only 
qualitative  agreement  is  seen.  This  is  due  to  the  non-realistic  constant 
eddy  diffusivity  model  used  in  the  numerical  calculations.  Numerical 
results  given  in  Figures  18  and  19  show  that  the  flow  and  temperature 
fields  are  more  constricted  about  the  centerline  for  larger  Gr^. 

In  reality,  the  turbulence  level  decreases  as  the  radius  and  the 
distance  from  the  plate  increases.  At  some  distance  outside  the  primary 
flow  region,  the  flow  field  (the  secondary  flow  due  to  entrainment)  is 
laminar,  and  the  description  of  the  flow  behavior  in  this  region  would 
be  determined  by  the  molecular  properties  (i.e.,  by  Gr).  There  must 
exist  some  sort  of  transition  of  the  governing  transport  mechanisms.  It 
appears  possible  to  improve  the  numerical  results  by  incorporating  a 
spatial  variation  of  Gr^  (i.e.,  the  eddy  diffusivity)  in  the  numerical 
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formulation  such  that  Gr^  would  increase  with  the  radius  and  the  eleva* 
tion  until  it  is  approached  to  the  ambient  value  Gr.  This,  however, 
would  require  an  extensive  numerical  experimentation.  Because  of  the 
large  amount  of  computer  time  which  would  be  required,  this  procedure 
was  not  pursued. 

Figure  44  shows  the  calculated  indraft  profile  for  the  case  #p  * 
520°F.  The  results  show  that  the  calculated  indraft  velocity  is 
consistently  larger  than  the  experimental  data.  The  percentage 
discrepancy  is  large  especially  at  large  elevations  (~90%).  Since  all 
measurements  of  indraft  velocity  were  in  the  laminar,  secondary  flow 
region,  it  is  conceivable  that  this  discrepancy  was  due  to  the  use  of 
too  large  a  value  of  Uc  for  nondimensionalization,  or  it  may  have  been 
due  to  the  inability  of  describing  the  whole  flow  field  by  means  of  a 
single  transport  mechanism. 


VIII .  CONCLUSIONS 


Based  on  the  findings  of  this  research  the  following  conclusions  are 
drawn  for  the  turbulent  free  convection  field  above  a  heated,  horizontal, 
circular  flat  plate: 

EXPERIMENTAL 

1.  The  turbulence  in  the  primary  flow  region  (Figure  4)  is 
characterized  by  low  frequency  (~10  cps)  and  large  amplitude  ( ~1  fps) 
turbulent  fluctuations.  The  fluctuating  velocities  are  the  largest  on 
the  centerline  and  decay  as  the  radius  increases. 

2.  Because  of  this  large  turbulence,  the  conventional  hot-wire 
technique  for  turbulence  measurements  on  linearized  theory  was  impossible 
to  apply. 

3.  There  exists  a  thermal  boundary  layer,  within  approximately 
0.02"  of  the  plate  surface,  in  which  the  major  part  of  the  overall 
temperature  drop  takes  place.  The  temperature  distribution  in  this 
region  may  be  calculated  based  on  molecular  conduction. 

4.  The  vertical  temperature  gradient  and  hence  the  local  heat 
flux  at  the  plate  surface  increases  significantly  with  the  radius.  Data 
show  that  the  temperature  gradient  at  plate  edge  is  approximately  twice 
the  value  at  centerline. 

5.  The  region  of  significant  deviation  from  ambient  temperature  and 
velocity  is  constricted  toward  the  centerline. 
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6.  The  induced  indraft  is  large  at  low  elevation  (~0.5  fps)  and  is 
largest  within  an  inch  of  the  ground.  This  indraft  is  on  the  sane  order 
of  magnitude  as  Uc  (a  characteristic  velocity  directly  obtainable  from 

a  dynamic  force  balance). 

7.  Hot-wire  method  has  been  developed  which  is  suitable  for  this 
type  of  velocity  and  turbulence  measurements  when  the  information  contained 
in  the  various  moments  of  the  hot-wire  voltage  output  is  used.  Details 

of  the  data  can  be  more  clearly  shown  as  the  order  of  the  moment  increases, 

being  a  more  sensitive  measure  of  the  phenomena.  However,  due  to  the 

amplification  and  accumulation  of  errors  during  the  process  of  talcing 

higher  moments,  the  proper  interpretation  of  the  high  moment  data  may  be 

a  difficult  task  (e.g.,  the  problem  of  separating  signal  voltage  and 

background  noise  may  be  an  analogy). 

S.  Numerical  simulation  of  the  the  hot-wire  voltage  output  is 

shown  to  be  a  useful  tool  for  the  actual  hot-wire  data  reduction. 

9.  The  eddy  diffusivity  can  be  estimated  either  from  the  energy 

2 

equation  using  basically  temperature  data  or  from  the  product  »ur  of  the 
turbulence  velocity  and  the  characteristic  time  data.  The  agreement  of 
available  data  is  within  15%  for  low  elevation  points  and  up  to  60%  for 
higher  elevations. 

NUMERICAL  CALCULATIONS 

1.  A  numerical  calculation  method  using  a  time  advancing,  explicit 
differencing  scheme  for  computing  essentially  the  laminar  free  convection 
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field  above  a  heated,  horizontal  disk  has  been  developed  and  validated. 
The  calculated  temperature  gradient  at  the  plate  surface  increases  with 
the  radius  (also  an  observed  experimental  fact).  The  numerical  heat 
transfer  results  obtained  agree  with  experimental  heat  transfer  correla¬ 
tion  (NuocRa*^)  to  within  4%  of  the  value  of  the  exponent. 

2.  The  centerline  temperature  data  were  compared  with  numerical 

results  and  the  best  fit  was  used  to  establish  the  appropriate  values 
of  Gr^  and  the  turbulent  Grashof  number  and  temperature  rise 

characterizing  the  flow.  The  match  between  calculation  and  experimental 
data  is  most  successful  when  Gr^  is  on  the  order  of  10*  (Grashof  number 
Gr  based  on  molecular  viscosity  is  on  the  order  of  10*0  for  this  work).  - 

3.  For  describing  the  gross  behavior  of  the  flow  field,  the  intense 
variation  of  eddy  diffusivity  in  the  boundary  layer  region  is  neglected 
by  the  specification  of  a  characteristic  temperature  rise  at  the 
plate  surface.  This  temperature  rise  is  determined  also  on  the  basis  of 
plate  heat  flux  and  eddy  diffusivity.  Only  qualitative  agreement  can  be 
obtained  between  the  numerical  and  the  experimental  data.  Clearly,  the 
constant  eddy  diffusivity  model  is  not  suitable  for  describing  the 
complete  flow  field. 


Figure  2.  Heater  installation  (approximately  half  actual  size) 
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Figure  6.  Thermocouple  inetallacion  in  the  heated  aluminum  place. 
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tunnel  -  schematic. 
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Figure  11.  Tunnel  velocity  calibrations. 
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Figure  12.  Hot-wire  calibrations  (0.001"  dia.  x  0.5"  long 

platinum  wire,  positioned  horizontally,  estimated 
wire  temperature  *  300°F) . 
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Figure  16.  Thermocouple  and  hot  wire  probe  arrangement. 
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Figure  18.  Horizontal  temperature  profiles  ($  «=  520°F) . 
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Radius,  Inches 

Figure  19.  Horizontal  temperature  profiles  ($  *  410°F). 
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Air  temperature,  F 

Figure  21.  Temperature  distribution  near  the  plate  ($  *  520°F). 
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Figure  23,  Comparison  of  wave  analyzer  outputs  for 

different  wire  lengths  (0  «  520°F,  X  ■  6" 
Y  ■  6",  wire  horizontal).  ^ 
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Figure  25.  Hot-wire  bridge  voltage  output  (Y  =6"  0  =  520°F1 
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Figure  26.  Autocorrelation  function  of  hot-wire  bridge  voltage  output  for 
centerline  locations  (wire  horizontal ,  6  •  520°F) 
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<d)  X  -  0,  Y  -  12" 
Figure  28.  (Continued) 
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Figure  29.  Indraft  profiles  (  X  *  12”  )„ 


Figure  32  (continued) 
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Figure  34.  Numerical  result*  at  various  stags*  of  calculations  for  a 
transformed  coordinate  system  (grid  system  *  17  x  17,  1-9', 
C1  *  C2  "  4»  6rT  "  7*43  *  l°9»  **  *  0.708). 


Figure  35.  Comparison  of  numerical  results  for  non-flow  condition  with 
exact  solution  of  diffusion  equation). 
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Figure  37.  Effects  of  domain  of  calculations  on  the  numerical  results  - 
vertical  extent  (GrT  =  10*.  Pr  =  0.708).  - 
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126 


X/D 

Figure  38.  Effects  of  domain  of  calculations  on  the  numrical  results  - 


horizontal  extant. 
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Figure  39.  Effect  of  mesh  spacing  on  centerline  temperature 
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Figure  41.  Effect  of  laesh  spacing 
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Figure  45.  Comparison  of  updraft  data  with  numerical 
results  (9p  *  520°F). 
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Appendix  A 

MOTIVATION  OF  RESEARCH 


The  origin  of  this  project  was  an  attempt  to  study  the  gross  flow 
behavior  :.n  the  vicinity  of  large  unconfined  fires.  This  problem  has 
been  the  subject  of  several  investigations  during  the  past  decade  (29, 
30,  31,  37).  Experimental  observations  suggest  that,  in  general,  the 

physical  behavior  of  a  large  fire  may  be  described  in  terms  of  the  four 
zones  shown  below  (56) : 


1. 


Plume  Zone 
Above  the  fire  th 
where  thia  curren 


ere  ia  a  rising  current  of  heated  gases.  The  zone 
c  is  the  dominant  feature  is  called  the  fire 
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convection  plume.  Entrainment  of  ambient  air  and,  therefore,  cooling 
of  the  rising  column  take  place.  The  flow  is  fully  turbulent.  Pire 
whirls  may  sometimes  be  present.  The  affected  arose  section  of  the 
surrounding  air  gradually  widens  along  the  path  of  the  ascending 
mass  of  the  hot  gases  until  its  momentum  is  balanced  by  entrained 
ambient  air.  The  flow  behavior  in  the  plume  zone  resembles  that  of 
a  free  jet.  The  boundary  layer  assumptions  are  valid  and  under  some 
circumstances  similarity  solutions  may  be  obtained.  Mathematical  and 
experimental  investigations  of  the  plume  have  been  fruitful  (e.g. , 
references  4,  5,  6,  7,  8,  52). 

2.  Secondary  Flow  Zone 

This  is  the  "upstream"  region  of  cold  air  outside  the  fire.  In  an 
otherwise  still  ambient  atmosphere,  the  fluid  motion  in  this  zone 
is  due  solely  to  the  pressure  defect  resulting  from  the  buoyancy  of 
the  hot  ascending  gases.  The  entrainment  and  the  associated  indraft 
are  largest  at  altitudes  relatively  near  the  base  of  the  fire.  This 
indraft  is  substantial  for  large  fires;  it  tends  to  prevent  the 
outward  lateral  propagation  of  the  fire.  Thus  when  the  indraft  is 
large  enough  all  the  burning  may  be  confined  in  a  relatively  stationary 
region.  This  is  one  reason  why  the  ground  indraft  is  so  important 
in  fire  behavior. 

3.  Surface  Interaction  Zone 

This  is  the  region  of  fuel  supplies  such  as  buildings,  trees,  etc. 

Major  destructive  burning  and  the  release  of  gaseous  fuels  take  place 
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la  this  zone.  Inside  this  zone  are  intricate  interactions  of  heat 
and  mass  transfer,  chemical  reactions,  and  air  entrainment.  In  the 
start-up  phase  of  a  mass  fire,  ignition,  spreading  and  the  interactions 
between  individual  small  fires  occur  in  this  zone.  The  burning 
process  in  this  region  is  usually  incomplete.  In  a  large  area  fire, 
this  zone,  although  the  source  of  energy  supply,  is  usually  confined 
in  a  region  very  close  to  the  ground  level  and  is  very  small  compared 
with  that  of  the  whole  region  which  is  significantly  influenced  by 
the  fire  (Reference  29).  It  is  believed  that  the  detailed  physico¬ 
chemical  phenomena  associated  with  combuation  in  this  region  are 
only  local  effects,  and  that  the  large  scale  convective  motion  in  a 
fire  is  governed  primarily  by  the  heat  release.  Thus,  for  modeling 
purposes,  we  assume  that  the  detailed  knowledge  of  phenomena  in  this 
zone  is  not  essential  on  understanding  of  the  gross  flow  behavior  of 
large  fires. 

4.  Core  Zone 

Near  the  base  of  the  fire  is  a  region  linking  the  surface  interaction 
zone,  the  secondary  flow  of  fresh  uir  and  the  convection  plume  above. 
Here  the  pressure  defect  driving  the  low  elevation  secondary  flow  is 
developed.  The  winds,  temperature,  and  oxygen  concentration  experi¬ 
enced  by  the  surface  interaction  zone  are  determined  by  the  core 
phenomena.  The  general  understanding  of  the  physics  in  this  zone  is 
essential  to  the  studying  of  fire  flow  fields.  The  flow  field  in  this 
zone  cannot  be  described  mathematically  by  the  boundary  layer  type 
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formulation  as  used  for  convection  plume  calculations  (4,5,6,7,8,52). 
This  aspect  is  demonstrated  by  Stewartson  (32)  with  an  example  of 
the  natural  convection  above  a  heated  horizontal  plate. 

The  fire  problems  studied  to  date  have  been  concerned  primarily  with 
the  convection  plume  well  above  the  fire.  Solutions  are  based  on  free 
boundary  layer  approximations  and  the  assumptions  of  small  density 
changes  and  similarity  profiles.  Results  of  these  calculations  are  in 
good  agreement  with  available  experimental  data  (4,5,7,52).  Experimental 
investigations  of  the  behavior  of  large  fires  have  been  carried  out  by 
personnel  of  Pacific  Southwest  Forest  and  Range  Experiment  Station. 

Wood  piles  or  simulated  houses  covering  areas  up  to  forty  acres  were 
burned.  Although  numerous  data  were  collected,  only  a  very  limited 
amount  of  fire  flow  field  data  have  so  far  been  made  available  to  the 
public  (29). 

Due  to  the  complexity  of  the  details  of  the  fire  problem  and  the 
prohibitive  cost  of  field  experiments,  it  is  desirable  to  utilize 
mathematical  models  to  the  greatest  extent  feasible.  A  logical  start 
is  to  try  for  an  approximate  description  of  the  large  scale  flow  features 
with  a  simple  mathematical  model,  and  to  then  attempt  to  fill  the  gaps 
by  adjustment  to  agree  with  field  data  as  available. 

As  discussed  earlier,  the  flow  field  generated  by  a  fire  is  driven 
primarily  by  free  convection  effects  resulting  from  the  heating  of  the 
air  in  the  proximity  of  the  fire.  We  regard  the  combustion  processes 
in  the  surface  interaction  zone  as  local  in  nature.  We  then  assume 
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further  that  the  heet  release  from  the  fire  is  the  dominant  factor  which 
govern  the  gross  behavior  of  the  fire  flow  field  so  that  the  fact  of 
combustion  may  be  modeled  by  prescribed  temperature  or  heat  flux  distri- 
bution  at  the  base  of  the  fire. 

A  limiting  case  of  this  is  the  free  convection  field  produced  by  a 
prescribed  uniform  heat  source.  Hopefully,  this  will  qualitatively 
describe  the  gross  features  of  the  flow  Induced  by  a  large  area  fire  and 
will  provide  a  starting  point  for  the  fluid  mechanical  study  of  a  fire- 
like  flow  field.  Once  the  problems  associated  with  this  fluid  mechanical 
model  are  overcome,  as  a  second  step,  a  more  realistic  model  including 
fuel  and  air  mixing,  followed  by  gas  phase  coeibustion,  could  be  devalopad, 
for  example,  by  simulated  spatial  distribution  of  sources  of  heat  and 
mass. 

For  mathematical  simplicity,  the  free  convection  due  to  a  horizontal 
heated,  circular  surface  flush  with  the  ground  surface  in  an  otherwise 
still  ambient  is  considered  as  the  model.  This  model  provides  a  flow 
field  which  describes  the  basic  features  of  a  large  fire  as  discussed 
earlier  in  this  appendix. 

The  free  convection  model  has  also  been  considered  as  an  approach 
to  the  scaling  of  the  indraft  velocity  due  to  an  area  fire  (31,37,57). 

By  a  dimensional  analysis.  Byram  (57)  has  shown  that  preserving  the 
Grashof  number  between  model  and  prototype  is  impossible  for  any  useful 
length  scaling  ratio.  Fortunately,  at  the  expected  magnitude  of  Gr  in 
a  fire  (~10A^),  the  flow  is  fully  turbulent  and  molecular  transport 
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phenosMna  have  been  postulated  unimportant.  Thus  If  Gr  can  be  neglected, 
for  geometrically  similar  models,  the  heat  flux  and  indraft  velocity  are 
proportional  to  the  square-root  of  the  length  scale.  The  validity  of  this 
concept  has  been  demonstrated  experimentally  by  Lee  (31)  as  conclusively 
as  possible  with  the  limited  data  available  (30). 


Appendix  B 
ENERGY  BALANCE 

The  calculations  for  the  energy  balance  terms  summarized  in  Table  1 
are  discussed  in  this  appendix.  Three  types  of  energy  loss  are  considered 
here: 

a.  Radiative  losses  from  the  aluminum  plate,  the  asbestos  ring, 
and  the  floor. 

b.  Convective  losses  from  the  asbestos  ring  and  the  floor. 

c.  Conductive  losses  through  the  insulation  jacket  and  the  floor. 

As  an  example,  the  calculations  for  Run  No.  1  are  described  below: 

1.  Radiative  Losses 

The  following  equation  was  used  for  all  radiative  loss  calculations: 

Qr  -  «rA(#4  -  Btu/hr 

where  #  (* 537°R),  *(■.1714  x  10’®  Btu/hr-f t2-R^) ,  «,  and  A  are, 
respectively,  the  ambient  temperature,  Stefan-Boltzmann  constant, 
emissivity,  and  heat  transfer  area. 

2.  Convective  Losses 

From  the  experimental  correlation  Nu  ■  0.141  Ra^®  (Reference  26) 
and  the  definition  of  qc  *  hAS,  we  estimate: 

qc  -  0.141k  (Pr  ^fj1/3A#4/3 
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Using  the  property  values  for  air  at  #^(58) ,  the  above  expression 
may  be  simplified  to: 

Qc  -  0.2434#  1,33  A  Btu/hr. 

The  calculated  energy  losses  for  the  above  two  items  are  summarized 
below: 


Part  Name 

Aluminum  Plate 

Asbestos  Ring 

Floor 

t 

.05 

.96 

.91 

A,  ft2 

3.18 

.84 

1.6 

6,  °R 

979 

740 

597 

Qr,  Btu/hr 

229 

294 

110 

Qc»  Btu/hr 

231 

86 

3.  Conductive  Losses 

These  losses  are  from  the  insulation  jacket  and  partly  from  the  floor. 
The  temperature  at  two  points  along  various  heat  transfer  paths  were 
measured  (see  Figure  2  for  thermocouple  locations).  The  total  heat 
transfer  surface  was  divided  into  several  parts.  For  each  of  these 
surfaces,  an  average  temperature  drop  was  idhtaiaed  from  eha  measured 
temperature  data.  Since  the  conduction  surface  area,  the  length  of 
path,  and  the  thermal  conductivity  (No.  4  vermiculite,  see  footnote 
on  p.  15  )  are  known,  the  heat  transfer  was  calculated  by: 

Q  •  -k  ~  A  Btu/hr. 

The  calculated  total  conductive  loss  through  the  insulation  jacket 
was  347  Btu/hr  and  the  conductive  loss  through  the  floor  was  45  Btu/hr. 


Appendix  C 

TEMPERATURE  CORRECTIONS 

The  temperature  indicated  by  the  thermocouple  is  corrected  for 
losses  due  to  radiation  and  conduction.  We  refer  to  the  detailed 
configuration  of  the  thermocouple  give  in  Figure  20. 

Heat  is  transferred  to  the  exposed  surface  (0.4"  long)  of  the 
thermocouple  wire  by  convection  and  lost  from  the  same  surface  by 
radiation.  Heat  is  also  lost  by  conduction  through  the  portion  of  the 
thermocouple  wire  in  the  stainless  steel  tubing.  For  thermal  equilibrium, 
we  require  that  qconv  =  qcond  +  qrad;  this  gives: 

Aconvh<*a  '  *t>  "  Acondkvlre  <*t  '  *  Ar.d  '*<4  '  '«> 

where 

Aconv  11  Arad  *  +  rD2  =  7.29  x  10'5  ft2 

Acond  "  Fd2/4  "  2-16  x  1C"7  ft2 

K^ire  *  ^copper  +  ^constantan  =  218  +  12.8  *  231  Btu/hr-ft-R 
««  emissivity  «  0.8 
<r«  .1714  x  10”8  Btu/hr-ft2-R4 
■  ambient  temperature  *  537°R 
9a  -  true  mean  temperature  of  the  air,  °R 
6t  *  mean  temperature  indicated  by  thermocouple  output,  °R 
L  «  length  of  conduction  path  -  0.25  ft  (see  Figure  20). 
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Substituting  the  appropriate  values  into  the  above  equation  and 
rearranging,  ve  get: 

r 

Lh 

noo 


a# 


,  .  2  75  + 

#t  h  T  h 


4 


(C-l) 


where  the  first  term  on  the  right  hand  side  of  the  equation  is  the 
correction  due  to  conduction  loss  in  °F  and  the  last  term  is  the  correction 
due  to  radiation  loss  in  °F. 

Consider  the  heat  transfer  surface  in  the  exposed  portion  of  the 
thermocouple  as  composed  of  a  sphere  of  diameter  D  (for  the  bead)  and  a 
circular  cylinder  of  diamter  d  (for  the  wire);  an  equivalent  wire 
diameter  weighted  by  surface  area  was  used  for  heat  transfer  coefficient 
calculations. 


4.  ■  qttVJpfr  •  «.«*«■. 


At  an  assumed  air  temperature  (e.g.,  that  calculated  based  on  conduction 
only),  the  properties  of  air  may  be  found  in  Appendix  III  of  Reference 
58),  with  the  local  mean  velocity  from  velocity  measurements,  the  local 
Reynolds  number  may  be  calculated.  From  heat  transfer  correlation  (e.g. , 
Reference  58,  Figure  7-3),  the  heat  transfer  coefficient  h  may  be 
calculated.  The  temperature  corrections  thus  obtained  are  shown  in 
Figures  21  and  22.  A  numerical  example  is  given  below. 

Consider  the  point  X  -  0,  Y  »  0.02"  for  Run  1  (Figure  21).  The 
local  mean  velocity  was  obtained  from  interpolation  between  the  points 
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Y  *  0  and  3".  Qa  was  estimated  to  be  458°F.  The  following  additional 
data  were  used: 

V  *  0.009  ft/sec 
#t  -  782°R 

k  -  .0224  Btu/hr-ft-P 
*  -  4.3  x  10"4  ft2/sec 
Re  ■  Vde/»  *  0.024 

since  this  Reynolds  number  is  so  small  (forced  convection  data  unavailable), 

the  actual  heat  transfer  is  in  the  free  convection  regime.  Setting  Ra  * 

2 

PrGr  ■  PrRe  and  using  Figure  7-3,  Reference  (58),  we  get: 

Nu  -  0.51 

h  -  Nu(k/de)  -  9.5  Btu/hr-ft2-F. 

Substituting  the  appropriate  values  into  Equation  (C-l),  we  get: 

M  -  116°F 

Finally , 

la  «  ^  +  Al  -  438°F 

Error  Analysis 

In  the  above  calculation  procedure,  we  need  to  estimate  the  values 
of  #a  and  V  for  computing  the  heat  transfer  coefficient  h.  From  the 
dependence  of  h  on  the  governing  parameters,  it  can  be  shown  that  h 
depends  only  weakly  on  the  value  of  ia  and  more  strongly  on  V.  Experi¬ 
mental  heat  transfer  correlation  (Reference  58)  shows  that  the  change 
of  Nu  is  very  gradual  in  the  range  of  our  interest  (Ra-vlO"4).  Spot 
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checking  of  the  results  were  made  and  it  was  found  that  an  error  of  20°F 
on  #a  would  result  1°F  and  a  20%  error  of  V  would  result  a  5°P  change  in 
the  temperature  corrections.  Since  20°F  is  larger  than  the  possible 
variation  in  our  temperature  measurements  and” 20%  error  in  V  is  a 
reasonable  upper  bound,  the  estimated  error  in  our  temperature  corrections 
would  be  less  than  6°F. 

The  source  of  error  during  data  taking  was  basically  due  to  the 
drifting  of  the  recorder  chart  paper,  <1  mm  is  a  quite  noticeable  amount). 
This  drift  would  represent  a  maximum  possible  error  of  0.1  mv  in  the 
thermocouple  output,  or  an  error  of  ~5°F.  Consequently,  the  total 
possible  error  in  our  temperature  measurements  would  be  approximately 
11°F.  Since  the  main  contribution  of  correction  is  due  to  radiation  loss 
which  increases  rapidly  with  the  wire  temperature,  it  is  clear  that  the 
amount  of  correction  and  hence  the  associated  error  will  be  greatly 
reduced  when  the  thermocouple  gets  away  from  the  plate. 


Appendix  0 

MUIVATKW  OF  EQUATION  (5-4) 


5F  -  <S2-A’)2  «  [<I2-A)  +  (A-A ' )j 2  -  (P-A)2  ♦  2(A-A‘)(P-A)+(A-A’)2 
Multiplying  by  CX  and  rearranging  the  firet  of  Equation  (3-4)  will  reeult. 
W'2  -  (E2-A’)4  -  |(E2-A)  +  (A-A 1 ) j 4 

-  (E2-A)4  +  4(A-A')(E2-A)3  +  6(A-A')2(P-A)^ 

+  4  (A-A')  3  (P-A)  +  (A-A’)4 

Since, 

(E2-A)3  -  [(E2-P)  +  (P-A)j 3 

-  <E2-E2)3  ♦  3(E2-A)(E2-E2)2  +  3(P-A)2(E2-P)  +  (P-A)3 

-  <E2-P)3  +  3  (P-A)  J(E2-A)  -  (P-A)]2  +  (P-A)3 

-  <E2-P)3  +  3<e2-A)  [(E2-A)2  -  (P-A)2]  +  (P-A)3 

-  (E2-E2)3  +  3 (P-A) (E2-A)2  -  2 (P-A) 3 

We  have, 

W'2  -  <P-A>*  ♦  4(A-A’)(E2-P)3  +  12 (A-A ' ) (?-A) (P-A)^-  8(A-A')(P-A)3 

+  6  (A-A 1 ) 2  (P-A)  ^  +  4(A-A')3(E2-A)  +  (A-A*)4.  (D-!) 

Multiplying  by  C^2  and  neglecting  the  tern  4Ci2(A-A')(E2-E^)3,  the  a^cond 
of  Equation  (5-4)  will  result. 
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This  neglected  term  is  not  zero  in  general  but  is  very  smell  relative 
to  the  retained  terms  for  an  assumed  probability  density  of  E.  To 
illustrate,  let  this  distribution  be  Gaussian, 

p(E)  -^7- exp  [-a2(E-E)2] 

o  2  #  i 

where  £  is  the  average  of  E;  a  ■  (2#£  )  and  »  variance  of  E. 


Letting  u  «  a(E-E)  we  have 

?  .£•  E2  PdE  -  ji  +  rj2  e-“2du  - 12  + 

(E2-?)3  -  f* (E2-E^)3PdE 

■  Jr  xr[i^ + ')2  •  ^ + 1?)] 3  e‘u2  du 

“  "JTrXo”  [u2  +  2)S  U  '  1/2]  3  *’u2  du 


-  js  <1  +  6^2  E2) 

or, 

(E2-?)3  -  8  e£4  (1  +  3E2) 

From  the  data  given  in  Table  6,  the  value  of  (E^-P)3  and  hence  the 
neglected  term  4Ci2C2(E2-E^)3  can  be  calculated.  Calculations  show  that 
this  term  is  on  the  order  of  10~3  -  10-3  and  its  adjacent  terms  in 
Equation  (D-l)  are  on  the  order  of  1.  Therefore,  the  neglect  of  this 
term  does  not  introduce  apprecialbe  error  and  is  well  justified. 


Appendix  £ 

RANDOM  VARIABLE  ALGORITHMS  - 
NUMERICAL  SIMULATION  OF  HOT-WIRE  DATA 

Two  basic  schemes  for  obtaining  the  random  variables  for  hot-wire 
data  simulation  are  described  in  this  appendix.  The  parameters  to  be 
specified  consist  of  the  standard  deviation  of  the  fluctuations  and 
correlation  coefficients  as  defined  in  Equations  (5-8)  and  (5-9).  In 
this  procedure,  a  set  of  four  random  numbers  are  generated  for  computing 
the  normalized  fluctuating  velocities  as  defined  in  Equation  (5-7).  The 
computation  involved  is  such  that  the  fluctuating  velocities  are  consistent 
with  the  specified  correlation  coefficients.  The  algorithms  for  these 
computations  are  described  below. 

For  convenience  of  discussion,  the  flow  variables  (7,  V,  \f,  T)  in 
Equation  (5-7)  are  redesignated  by  4^,  i  *  1,  2,  5,  4;  thus  the  correla¬ 
tion  coefficients  become  4j_4j.  The  problem  is  to  compute  a  random 
sequence  of  sets  of  4^  consistent  with  the  specified  set  of  values  of 
♦i*j-  Two  schemes  for  obtaining  4^  are  described  below. 

NORMAL  DISTRIBUTION  SCHEME 

Let  4ji  be  independent  random  variables  with  zero  mean  and  unity 
standard  deviation.  Let  F  be  the  probability  distribution  function  of 
the  flow  variables  4^,  and  set: 
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An  approach  Co  smasurement  problems  based  on  principle  of  luit 
bias  has  baan  discussed  by  Leas  (55)  and  Janas  (59),  Poi lowing  the 
formalism  presented  by  Lees,  we  find  chat  if  only  first  and  second 
moments  of  the  flow  variables  (u,  v,  w,  #')  are  prescribed,  e.g . ,  \T,  v, 
»u2,  uv,  etc.,  the  least  biased  distribution  function  consistent  with 
the  given  set  of  moments  and  correlation  coefficients  is  given  by 
Equation  (E-l)  with  the  ♦' s  Gaussian.  Gaussian  distribution  was  chosen 
for  all  basic  calculations. 


AN  ALTERNATIVE  GENERAL  SCHEME 

This  scheme  is  the  same  as  above  except  for  higher  correlations  and 
is  designed  for  consideration  of  a  number  of  random  number  probability 
distributions.  For  a  given  set  of  correlation  coefficients 
the  following  set  of  relations  are  first  written: 


*2*3  "  bl*L*2  +  b2 

Mi  "  bl  +  b2  K*2  (E-4) 

♦l%  "  ci  +  c2*1*2  +  C3Vfj 
*2*U  m  Cl*fc*l  *  C2  +  Cl*2*l 
*1%,  “  Cl*3*l  +  c2*l*2  +  c3 
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where  a,  b^,  b2,  c^,  c2,  c3  are  constants  Co  b«  solved  la  tarai  of 
The  confutation  will  proceed  If  these  coo a cant a  satisfy  the 


follow lag  requirements: 


|.|  SI 

Ibj]  ♦  jb2|  S  1  (B-5) 

|cl|  *  c2|  +  |c3|  *  1 


Otherwise,  the  set  will  be  celled  inconsistent  and  the  confutation 
is  ce ruinated . 

Take  four  independent  randan  nunbers  ^,1*1,  1,  3,  4  fron  a  set 
of  specified  probability  density  (e.g,  Gaussian  distribution)  and  four 
swre  independent  randan  nunbers  R^,  i  *  1,  2,  3,  4  fron  a  set  which  is 
uniformly  distributed  on  the  interval  <0,  1)  and  define  G(R,  c)  such 
that: 


G(R,  c)  -  1  for  R  S  c 
-  0  for  R  >  c 


(E-6) 


The  flaw  variables  ^  are  then  calculated  as  follows: 


*2  -  5ign(a)G(R1,  Jaj)^  +  G(R]L,  l-|a|)*2 

♦3  -  Sign(b1)G<R2,  jbj)^  +  Sign(b2)G(R2,  jb2|)d2 

+  c<r2.  l-K|-h|>*3 


(E-7) 
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\  -  Sign(c1)G(R3,|c1|)^1  +  Sign(c2)G(R3, fc2|)*2  +  Sign(c3)G(R y 
|c3|)*3  +  G<R3,  1-|cJ-1c2|.Jc3|)#4  (E-7) 

GENERATION  OF  THE  »j 

The  random  variables  ^  used  above  are  generated  from  two  basic 
sets:  a  normally  distributed  set  with  mean  zero  standard  deviation  unity 
and  a  uniformly  distributed  set  R^  on  the  interval  (0,  1).  The  former 
was  used  for  obtaining  the  normally  distributed  ^  and  the  latter  was 
used  for  obtaining  the  normally  distributed  R^  and  three  other  types  of 
distribution  of 


#  -  2  VT* (R  -  0.5) 

♦  “  +1  for  R  £  0.5 

-1  for  R  <  0.5 


(E-8) 


(E-9) 


(E-10) 


Equations  (E-8)  and  (E-9)  describe  a  uniform  distribution  of  t  in  the 
range  -  \A  S  ♦  and  an  equal  distribution  of  #«*akl,  respectively. 
Equation  (E-10)  is  derived  from  Cauchy  distribution. 


/(♦)  *  |  (1  +  *2)"2 

r.LjL£UJL. 

I  /(♦')  d  * 
•'-00 


Let 
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whence 

*  (R  -  1/2)  *  +  tan  ^  ,  |tan  J 

Expand  the  above  equation  for  *»0  and  rearrange, 

*  (1  -  R)  -  3^  (‘A  -  sfy  ) 

Solving  this  equation  for  *,  Equation  (E-10)  ia  obtained.  Equation 
(E-10)  ia  uaed  for  the  tail  end  values  of  R  (i.e.,  for  R»Q).  A 
rejection  scheme  baaed  on  Equation  (E-10)  ia  described  in  the  following 

block  diagram: 
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COMPUTER  PROGRAM 

This  program  is  written  in  FORTRAN  IV  computer  language.  Three 
permutation  schemes  have  been  used  in  the  main  program  for  the  purpose 
of  parametric  studies.  Since  there  is  no  difficulty  in  programming  these 
schemes,  only  one  of  them  is  included  here. 

Input  Data 

1.  AW,  BW,  CC,  CN,  AL,  GAM  correspond  to  A,  B,  k2,  0.5,  «,  r,  in 
Equation  (5-10). 

PHIV,  THEV,  VM  correspond  to  *v,  9 v>  Vm  in  the  sketch  on  p.  45 
ID  -  display  index:  >2  means  complete  display;  -1  means  standard 
display. 

NO  -  a  number  to  initialize  random  numbers. 

2.  SUTI(I) ,  SVTI(I) ,  SWTI(I)  correspond  to  a  matrix  of  values  of  »u, 

V  <v 

3.  CKUV,  CKVW,  CKWU,  CKUT,  CKVT,  CKWT  correspond  to  the  six  correlation 
coefficients  in  Equation  (5-9). 

STT  = 

DPHIW  »  A*w  (interval  of  *w>  wire  orientation) 

Output  Data 

1.  UBX,  VBX,  WBX,  TBX  are,  respectively,  the  simulated  average  values 
of  u,  v,  w,  . 

2.  SUTX,  SVTX,  SWTX,  STTX  are,  respectively,  the  simulated  values  of 


V  V  V  v 
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3.  CKUVX,  CKVKX,  CKWUX,  CKUTX,  CKVTX,  CKWTX  are,  respectively,  the 

simulated  values  of  K  .  K  ,K  ,K.K.K. 

uv  *  w  wu’  u#  v#  ’  w# 

4.  EL,  LOB,  LIB,  L2B  are,  respe.  tively,  the  simulated  values  of  E, 
wl/2,"y,  w2#  gE,  SLO,  sLl,  SL2  are,  respectively,  the  standard 
deviation  of  the  simulated  values  of  E,  W 1/2,  W,  W2. 

5.  RLO,  RL1,  RL2  are,  respectively,  the  values  of  wl75,  w,  W^"  divided 
by  their  corresponding  values  for  =  0°. 

6.  R01  *  (Wl/2)2/w,  R02  =(wl/2)4/w2 ,  R12  =  W2/W2. 


HOT  WIRE  SIMULATION  PROGRAM  -  PARAMETER  VARIATION  SCHEME  NUMBER  1A 
SYSTEMATIC  VARIATION  OF  SINGLE  VARIABLE  DISTRIBUTION  FUNCTIONS 
E*SORT( ( ) .0+AL*TT)*AW+(l • 0+GAM*AL*T  T ) *RW* ( SORT i UPP2+CC*LPL2 ) )**CN I 
REAL  L0B.L1B,L2B.LCR.L1R.L2R 
LOGICAL  HAG 

0 1  MENS  1 01.  Ct4,4).B(4»4,4),SC<4>»S«(4.4)»CK(4»4l .SUT  l ( 100) *SVTI ( 100 
1I.SWTI ( 1  GO ) »  SUMF ( 20 ) »UTT ( 5000 ) »VT  T ( 5000 ) » WTT ( 5000 ) »TTTI 5000 ) »Z (4 ) * 
?R(3) 

COMMON  AW.BW  »AL  »CC  *CN .PHI V. THEV.PH 1 W» THEW. VM.SUT  »SVT*SWT*STT  »CKUV  * 
1  CK  VW « CK Wl  »CKUT  .CKVT .CKWT  «NmAX .IC.ID.CK.R.SB.C.SC.FlAG.EB.SE .LOB.SL 
20.L1B.5L1  »L?B  »  SL?  *OPH  I  W.Z  »R«  SljME.  I  IT  T  »VTT  »WTT*  TTT  *GAM 
REAO (5*1  )AW.HW,CC.CN,AL*GAM,PriIV.THEV.VM,lD,NO 

1  FORMAT <2K8.4.4F7.4.2F8.3*F8.4/I2»I6I 

WRITE '6.?)AWtRW,CC.CN»AL.GAM,PHlV,THEV.VM,ID,N0 

2  FORMAT ( 1  * 3H1H0T  WIRE  SIMULATION  PROGRAM  -  PARAMETER  VARIATION  SCHF 
IMS  1A-  VARIATION  OF  OIST.  FCN.  INDEX  IC  AND  NO.  TRIALS  N  /67H  E«S 
20RT ( ( 1 . 0+ AL«  TT ) *AW  + (1 . 0+GAM«AL  *T  T 1  *  BW* ( SORT ( UPP2+CC*UPL2 ) ) **CN ) /4H 
3GAW=F7.4.2X3Hfc!W=F7.4»2X3  :C*F7 .4, ?X3HCN*F7.4.2X3HAL«F7.4.2X4HGAM*F 
47.4.2X5HPHlV*F8.3,?X5KTHEV*F8. 3.2X3HVM*F7.4/4H  I D» I  2. 2X3HN0* 1 6 ) 

RN*RNU ( NO ) 
rm.RHNINOI 
1*0 

3  1*1  +  1 

READ ( 5 .4  I SUT I ( I ) .SVT I ( 1 ) *SWT1  (I) 

4  FORMAT ( 6F 1 0 « 5 ) 

IFISUTI ( I I.GE.O.OIGO  TO  3 
!$*!-l 

5  READ! 5.4 iCKUV.CKVW.CKWU.CtCUT.CKVT.CKWT.STT.DPHtW 
IFICKUV.LT.O.OlGO  TO  10 

WRITEI6.fc )CKUV,CKVW,eiCWU.CKUT,CKVr.CKWT.STT,DPHlW 

6  FORMAT ( 6H0CKUV=F7.4  *3X5HCKVW=F7»4  »  3X5HCKWU*F7 .4.3X5HCKUT *F7»4.3X5H 
1CKVT*F7.4.3X5HCKWT*F7.4.3X4HSTT*F7.3,3X6HDPHIW«F8.3) 

DO  9  1*1,  IS 
SUT  *  SUT I  I  I  ) 

SVT *SVT I (  I  ! 

SWT* SWT  I (  I  i 
WRITE<6«7)SUT,SVT,SWT 

7  FORMAT(5HOSUT*F9.5.5X4HSVT*F9.5.5X4HSWT*F9.5) 

DO  9  J*  1.4 

IC*J-1 

DO  9  K*1  ,3 

I F ( K • EQ« 1 ) NMAX* 100 

IF(K.FO.L1NMAX=400 

IF(K.EQ.3)NMAX*1000 

WR I TE ( 6  » t ) IC.NMAX 

8  FORMAT (4hOIC«I?,5X5HNMAX*I4) 

CALL  XWIPF 

I F ( . NOT . FLAG ) GO  TO  5 

9  CONTINUE 
GO  TO  5 

10  STOP 
END 
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SUBROUT  INF  XWIRF 

REAL  L0B.l18.L2B.L0R.L1R.L2R 

logical  Flag 

DIMENSION  CU.A  )  ♦B(fc.A.A)  •  SC «  A)  .  SR (4,4)  .CKI4.4)  .SUTI <  100) .SVTI <  100 
1 ) ,  SWT  I ( lf.OI .Sl)ME(2C).UTT(SOOC) .VTTtSOOO) «WTT ( SOOC ) .T1 TI 5000 ) ,Z (4 > « 
79(31 

COMMON  AW.8W.AL.CC.CN.PHIV.THEV.PhIW.THEW.VM.SUT  .  SVT . SWT »ST T .CKUV ♦ 
ICKVW.CKWU.CKUT.CKVT.CKWT.NMAX. IC.1D.CK.R.SB.C.SC.FLAC.EE.SE.L0B.SL 
20.L18.SL1.L2B. SL2.DPHIW.Z.R. SuME.uT  T *VT  T .WTT.TTT .GAM 
THEW»0.0 
PHlW«-DPHlW 

1  PHIW«PH1 W+0PH1W 
IFfOHIW.GT. 180. 0100  TO  2 

IF !  PHIW.GT .90 .0. AND.CKUV.EO.O.C. AND.CK VW.EQ.O.O. AND.CKVT .EQ.P.O.AN 
ID. CXWT.EO. 0.0. AND. CKWU.EO. 0.0)00  TO  2 
CALL  WI9FSP 
IFI.NOT.FlAGIOO  TO  3 
GO  TO  1 

2  IF<CKUV.EO.CKWU.AND.CSCVT.EO.CKWT)GO  TO  3 
THEW-90.1 

PHIWO.O 
CALL  WIRtSP 

3  RETURN 
END 
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SURR0UT1N"  VJPFSP 

RFAL  LOB.LlB.L2B.LOR.LlR.L2R 

logical  flag 

DIMENSION  C<4,4  )  .8!4.4,4).SC<4).SN(4,4J  .CICI^.A!  .SUT1  mOC  I.SVTIUOO 
1 ) .SWTI < 100! .SUMEI20) .UTT(5000| .VTTiSOOOJ .WTTISOOO) .TTTI 5000 > »Z 1 4) • 
2R  (  3  )  *  ALF  (4.4)  •  SUMA  (4.4) 

COMMON  AW,8W.AL.CC.CN.PHIV.tHEV.PhU'.THEW.'  M.SijT  »SVT,$wT*STT .CKUV. 
KKVtX.CKWU.CKUT.CKVT.CKwT.NMAX.  1C.  !D.CK.B,SB.C.SC.FLAG,£8.SE.L0B.SL 
20,Li8.SL 1  ,u2ft.SL2.DPHlw.Z.R.SUME.isTT  ,VT T ,*TT , TTT . ALF .  SUE  A*6AM 
IF (NMAX.LT.O (GO  TO  10 
FLAG-. TRUE. 

C<  (  1  «2)«CKUV 
ClC(2.3)«CKVW 
OC(  1.31-CKWU 
CK I  1 *4 ) -CKUT 
CK ( 2.4 ) -t KVT 
CK ( 3 .4  )  ■ CK  WT 
DO  1  1-1.4 
DO  1  J-1,4 

I F  (  I «  EQ<  J  >CK I  I .jl-1.0 
IF (  I  .GT.  I ; CK ( I ,J)-CK( J.I  ) 

1  I F  (  ABS ( Ct.l  I  * J I  )  .GT.  1,0) FLAG*. FAL$f« 
do  u  x«) .4 

scm-r.r 

DO  3  1*1 .4 
C (  I  ,K)«0.0 
SB  (  I  »K  )  »t  .0 

DO  ?  J«1 

B< I  ,J,K  )  >0.0 

IF( I .CQ. J.OP.J.FO.X.OR.X.FO.I )GO  TO  2 

8! 1 »J.K I , (CK l I . J)-C< ( J»K >*CKI X • I ) )/ ( 1 ,0-CK ( J.K )**2 ) 

B< I .K.JI -ICK I  I ,K l-CK(K,J)#CK< j.I ) 1  '(  1  •  0“CK  I  K  *  J  )  **2  ) 

IF  (  ABSIB  I  I  .J.KD+AHSI  A(  I  »K.  JM .GT. 1.0) flag*. FALSE. 

L-10-I-J-K 

C  I  !  .K  )  *  <  CK  (  [  ,K  )  *  (  1  ,."-CK  <  L»  J  )*»2  )-CK  ( L  »< I  * (CK (  I  *L  )-CK  I  L.J)*CK( J.I  )  ) 
1-CK (J<K)*(CK(J. |  !-CK(  I  »L)*CK(L»J>  I  l/i 1 .0-CK ( L . J ) *»2~CK ( I .L ) * ( CK I  I . 
2L)-CK(L.J)«CK(  J.I  l)-CK<  J.n*(CK(J,I  l-CKI  I.LIKKIL.JI)  I 

2  CONTINUE 

3  SC (K )»SC ( K)*ABS(C (  I  »K  )  ) 

if<  sc i< i ,gt, i,o iflag -.false. 

4  OONTIN'JF 

IF! ID.GT.-2.0R.FLAGIG0  TO  6 

WRITF(6.S)i(CKlI.J),J«l»4).I«1.4!,(((B(I.J.K),K*l.A).J*1.4),I-1.4) 
1 • <  <  C( I ,J)  »J«  1  .4  I  , 1-1.4), i SC  I  I  ), I-),4) 
b  FORMAT ( 1h0« ( 4F]2.4,SX,4Fi2.4)  > 

6  IF (FLAG) GO  TO  A 
WRI Tf (6,7) 

7  FORMAT! 37HOCORRELATION  DATA  INCONSISTENT  -  EXIT) 

RFTuRN 

8  IF! IC.GF.DIGO  TO  80 
DO  81  1-1,4 

DO  81  J- 1,4 
81  Al  F ( I , j ( >0, 0 
DO  83  N-  !  »?o 
DO  82  I-..4 
DO  82  J-l  ,4 
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SUMA< I *J)»0.0 
do  *,?  k-i.a 

82  $UMA<  t  «J  )  -SUMAI  1  .  J » ♦  At_F  <  1  r<  )*ALF  (  <  ) 

00  83  1*1. A 

00  83  J*  1  * A 

ALP! I *J> *0.5*1 -CK< I » J ) +  SUMA ( ItJl) 

IF ( I .£0.J)ALF< 1 .Jl-ALF ( 1 .J) +0.6 

83  IF  UD?EQ.-2)  WRITE!  6.84) ! ( AlF( I »J) »J*1»4I .1*1.4) 

84  FORMAT ( 9H0ALF ( l * J ) / 1H  .16F8.5) 

go  to  io 

80  A*CK  ( 1*2) 

61-CK (2.3) 

87*0.0 

IF<CK<l*2i**?.LT.l»0)Rl«A<3»l»?) 

IF< CK (1*2) **2.LT.1.0)B2»B <3.2.11 
n«CK<1.4) 

C2*0.0 

C3*0«0 

IF ( CK ( 1  *  2 ) *Ck  •  2  «3)*CK  <3»1>.NE.1.0)C1*C(1*4> 

IF(CK< 1*2 )*CK(2.3)*CK < 3.1 • »NE. 1.0 )C2*C { 2*41 

IF(CK(1»2)*CK(2.3)*CKI3»1).NE.1.0)C3*C(3»4) 

IF< lO.EO.-2.ANO.NMAX.GT.fi 1WRI Tfl 6.9 ) A. B1 .B2.C1 .C2.C3 
9  FORMAT ( 3H  AIF7.4.3X3HB1*F7.4»3X3HR2*F7 «4*3X3HC2*F7»4* 3X3HC2*F7»4  * 3 
1 X3HC3“F7  « A ) 

10  BB-PHIW/57.2958 
BETA-PHIV/57.2958 

ALA«<THE\-THFW)/57.2958  . 

E1-SIN(BE»-BETA)+C0S(BB)*SIN<BETA)*( l.O-COSI ALA) ) 

E2*SIN<ALA)*C0S(BRI 

E3-C0S< BB-BETA) -C0S< B81*C0S (BETA )*<1.0-C0S< ALAI) 

N*0 

DO  11  1-1.6 

11  SOME < I) *0.0 

IF<PHtW.UF,O.O.OP.THFW.NE.O.O)GO  TO  1A 
DO  12  1-7.20 

12  SUMF<n-0.0 
IFINMAX.l  T.mr-O  TO  14 

13  FORMAT! 1H0.5X2HUT.10X2HVT.10X2HWT .10X2HTT .10X2HE 1 . 10X2HE2. 10X2HE3 • 
1 9X4HUPP2 .8XAHUPL2.9X1HF) 

14  r>r\  ?fl  iNxl.lon 
N-N+l 

IFINMAX.LT. 0)GO  TO  17 
DO  16  1-1.4 
IF<  IC.LE.01ZI I)*RNN»01 
IF<  IC.GT.O>Z< I1*RNU<0I 
IF< IC.EQ.-l 1G0  TO  16 
IP<I.LE.3)R(T )*RNU(0 ) 

IF( IC.EO.l )Z< I ) -3 . A6A 1* ( 2 ( I) -0.5) 

I F  <  I  C.E0.2 )Z  < I »*5I6N(1.0.<2< I  1-0,3) ) 

IF <  IC.NE.31G0  TO  16 
ZA-AMIN1  <ZU  1  .  <1.0-211)1) 

IF < ZA. LT. 0. 02027)2  U 1-SJGN < (SORT ( 0, 355773*2 A** ( -2. 0/3.0 1 -0 • 79 ) 1 . (Z 


1111-0.5)1 

IFfZA.LT. 0.02027) GO  TO  16 
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15  RZ-RNUiO) 

Z(I)*4.0*(RZ-0.5) 

RZ*RNU<  0 ) 

!F<RZ.GT.il.0fZ(I)**2J**(-2>>G0  TO  15 

16  CONTINUE 
JFdC.GE.OIGO  TO  160 

UT»(1.0-„LF(  1,1  )  )*Zd>-ALF(l,2>*Z(2)-ALF(l,3)*Z<3)-ALF(l»4)*Z<4) 
VT«(1.0-mLF(2,2) )*Z<2>-AIF(2,3)*Z(3)-ALF(2,4>*Z(4)-ACF{2»1 1*1111 
WT*( 1»0-ALF ( 3*3 1 )*Z  1 3 ) -ALF ( 3*4  )*Z (4 ) -ALF { 3* 1 ) *Z  < 1 1 -ALFI 3  »2  >*Z  <  2  S 
TT*(1.0-ALF(4,4) )»Z(41-ALF( 4*1)*Z(1 > -ALF ( 4,2 > *Z <2 1 -ALF( 4»3 >*Z (3) 

GO  TO  161 

160  UT«Z (  1  ) 

VT*Z ( 2 ) 

1P(R<1)  .LP.ARS(A) )VT«UT*$rGN(l *0,A) 

WT«Z ( 3 ) 

IF|R(2).LF.ABS(B1 ) )WT«UT*S1GN(1.0,B1) 

!F(R(2)*LE«ABS(B1 1 *ABS ( B2 ) • AND*R( ?'.GT.ABS(B1 ) )WT«VT*SIGN( 1.0*B2) 
TT*Z ( 4 1 

IF(R(3) *LF.ABS(C1 ) )TT»UT«SIGN( 1.0, Cl ) 

IF(R<3) *LE.ABS(C1 )+ABS(C?).AND.R( 3).GT.ABS(C1 ) )TT»VT*$IGN( 1.0,C2) 
IF(Rt3).LF.ABS(Cl ) +ABS « C2 >♦ ABS <C3 ) . AND.R ( 3  I *GT . ABS( Cl  M-ABS 1 C2 1 ) TT» 
1WT*SIGN( 1.0.C3) 

161  UTT ( N 1 »UT 
VTT ( N ) *VT 
WTT ( N ) »WT 
TTT(N)«TT 
GO  TO  18 

17  UT»UTT(N) 

VT«VTT(N) 

WT  «WTT ( N 1 
TT-TTT(N) 

18  UT «IJT*5UT 
VT«VT*SVT 
WT  «WT  *SWT 
TT  »TT*STT 

IF(PHIW.NE,.0.0.OR.THEW.NE»0.0)G0  TO  19 

SUME(7)*SUME(7)«-UT 

SUME (81*  SUME ( 8 1 +VT 

SUME  (91®  SOME  (  9  )  +-WT  , 

SUMEI iO)rSUME( 101+TT 
SUME (11) sSUMF (11) +UT  *VT 
SUME (12) «SUME (12) +VT  *WT 
SUME (13) “SUME (13) +WT*UT 
SUME (  14 ) «SUME ( 1 4 ) +UT  *TT 
SUME  (15)  -»5UME  (15)  ♦VT*TT 
SUME  (  16  )  *-SUME  (  16)+WT*TT 
SUME  (17)  t-SUME  (17)*UT**2 
SUME( 18 ) »SUME ( 18 )+VT**2 
SUME ( 19 ) «SUME ( 19)*WT**2 
SUME  ( 20 )  *SUME  (20)+TT**2 

19  UPL2* ( (Vf+UT)*E1+VT«E2*WT*E3)**2 
UPP2«< VM*UT)**?+VT**2+WT**2-UPL2 

ESQ* ( 1 .0+AL*TT ) *AW+ ( 1 .0+GAM*AL*T T  >*SW* ( SORT  < UPP2+CC*UPL2 ) 1 **CN 
F.0.0 

IF(ESO.GT.O.O)E=SORT(ESO) 

CL*IESQ-AWI/BW 
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SUMF< ! >«SUMFt ] ) +F 

SUMP<2)*SllMF(?|+FS0  "■».  . 

SUMS  <  3 !  •  SUM£  {  3  )  ♦CL 
SUMF ( 4 ) *  SHMF ( 4 ) L  *•? 

SUM£  (  3  )  *  SI  IMF  i  3  )  +CL  **4 
Sump  (  m  «  si *mf  (  a  t  i  #*fl 

20  I F ( ID.EQ.-2.AND.NMAX.GT.0>WRITF<6,21 >UT »VT .WT *TT *E1 .E2*E3«UPP2*UPL 
l?*e 

21  FORMAT ( 1H  10E12.S) 

IF(N.LT.IABSINMAX)  >  GO  TO  14 
FN*FLOAT<NI 

NMAX»- I ARS ( NMAX ) 

F4»SUM?Ji)/FN 

E2B*SUME(2)/FN 

LOB«SUMF(3)/FN 

L1B«SUME(4)/FN 

L?B«SMMF(S»/Fw 

IF(PHIW.MF.O,O.OR.THEW.NE.O.O)GO  TO  24 

tj"X*SUMP(7)/PM 

VAX*  S'  IMF  l  A  J  /F*i 

VAX*SUME(9)/FN 

TBX*SUMF(10)/FN 

UVBX*SUMfm  )/FN 

VWBX*SUMr ( 1 2 1 /FN 

WUBX*SUME<13)/FN 

UTRX*SIJME(14)/FN 

VTBX*SUME(15)/FN 

WTBX*SUME ( 16 ) /FN 

U2BX»SUME(17)/FN 

V7AX*SUMFI 1A) /FN 

W2BX*SUME(1<J)/FN 

T2BX*SUME  t  20 ) /FN 

SIITX.-1  .A 

IF(U2BX.GT.UBX**2)SUTX«SQRT(U2BX-UBX**2> 

SVTX  — 1.0 

IF(V2BX.GT.VBX**2)SVTX«S0RT(V2BX-V8X*»2) 

SVTX--1 .0 

!F{W2BX.GT.WBX**2) SWTX*50RT (W2BX-wBX**2) 

STTX»-1.0 

IF{T2AX.GT.TBX**2ISTTX*SORT(T2BX-tBX**2) 

CKUVX* ( UVBX-IJRX*V8X ) / ( SUTX*SVTX ! 

C<VWX* ( VWBX-VRX»WBX I / <  SVTX*SWTX ) 

CKWUX* ( #UBX-WBX*UBX ) / ( SWTX*SUTX ) 

CKUTX*(UTBX-ljqx*TRX  I  / ( SUTX^STTX S 
CKVTX* ( VTBX~VBX*TBX l/(SVTX*STTX! 

CKWTX*(WTBX-WBX*TnXI/(SWTX*STTX) 

WRITE<6»22)UBX.VBX*WBX»T8X*  SUTXtSvTX  »$WTX* ST TX  tCKUVX  tCKVWX • CKWUX  «C 
lxuTx.rxvTx*rxwTx 

22  FORMAT ( SH0UBX*F9.4 ♦ 3X4HV8X*F9«  4*  3x4HW8X*F9.4«3X4HTBX*F9*4t  3X5H5UTX 
1«F8.4,3X5hSVTX*F8.4.3X5H$WTX*F8.4,3X5H$TTX*F8.4/7H  CKUVX-F7.4.3X6H 
2CKVWX*F7  >4* 3X6HCKWUX*F7.4« 3X6HCKUTX* F7.4 « 3X6HCKVTX*F7.4«3X6HCKWTX* 
3F7.4) 

1F(PHIW.NF.O.O.OR.THEM.NE.O.O)GO  to  24 
WRI TE ( 6  »  23 ) 

23  FORMAT ( 1H0 *3X4HPHI W »3X4HTHEW» 7X2HFB  *  5X2HSE  *6X3HL0R*4X3H5L0*6X3HL1B 
1*4X3HSL1*6X3HL?B*4X3HSL2»1X7M LOB/LOR  »1 X7HL 1B/LIR  * 1X7HL2B/L2R»5X3HR 
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101 .5X1HR02.5K5HR121 
lor«lor 
LIR-Llfl 
L2R-L2B 

24  SE*-1 .0 

IF(E2S.GT.EB**?tSE«SQRT< < E2B-EB**2 > /FN > /ABS( E8 > 

SL0«-1 • 0 

IF(L1B.GT.LOB**2ISI.O*SORT( <L1B-LOR**2>/FN>/ABS«LOB) 

Sll*-1.0 

IF(L2B.GT.L1B**2)SU»S0RT<  U28-L1B**2»/FN>/ABS<L1B» 

SL2«-1.0 

IF(SUME(t)/FN.GT.L?B**2»SL2«SORTf (SUME(6)/FN-L2B»»2)/FN)/ABS(L2Bt 

RLO-LOB/LOR 

RL1«L1B/L1R 

RL2*L28/L2R 

R01»L0B**2/L1B 

R02«L0B»»4/L2B 

R12«L1B**?/L2B 

V'RITE(6*25)PHIW*THEW»EBtSE*L0B.SL0#LlB*SLl«L2B*SL2tRL0*RLltRL2*R01 
1 *R02*R12 

25  FORMAT ( 1H  .2F7.2.4<F9.4*F7.4},6F8.4) 

RETURN 

ENO 


Appendix  F 

NUMERICAL  CALCULATION  PROCEDURE  AND  COMPUTING  TIME 


The  overall  numerical  calculation  procedure  is  outlined  in  the 
following  block  diagram  with  a  step  by  step  description: 
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The  details  of  the  above  block  diagram  are  described  below: 

1.  The  basic  input  data  consist  of  the  grid  system,  source  size, 

Gr,j.,  Pr,  convergence  criteria,  instructions  for  decision  making,  etc. 

A  list  with  explanations  for  these  input  data  is  given  in  Appendix  G. 

2.  Several  decisions  are  made  in  this  step: 

(a)  For  a  starting  calculation:..  Go  directly  to  step  3. 

(b)  For  continuation  of  calculations:  Read  input  data  cards  .< 

I 

from  intermediate  results  of  previous  calcula Lions.  ' 

(c)  For  subdividing  mesh  spacing:  Read  input  data  cards  from 
results  of  the  converged  solution  of  the  previous  calcula¬ 
tions;  and  set  the  approximate  values  at  the  added  mesh 
points  equal  to  the  arithematic  average  of  the  known  values 
of  the  neighboring  points. 

3.  At  t  =  0,  set  all  dependent  variables  equal  to  zero  except 
T^  =  1  for  x  ^  D/2. 

4.  Calculate  all  the  repeated  constants  in  the  computer  program. 

5.  Calculate  At  from  the  stability  criterion  (Equation  6-10). 

6.  If  the  execution  time  limit  is  exceeded,  go  to  step  15  for  card 
output.  This  card  output  will  be  used  for  future  continuation 
of  calculations  as  described  in  step  (2b). 

7.  Calculate  from  Equation  (6-11)  for  all  (i,j)  interior 

to  the  boundaries  (hereafter  these  points  will  be  referred  to  as  the 
interior  points  as  opposed  to  the  boundary  points)  and  then  calculate 
T1  ^  j  for  all  interior  points  as  follows: 
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T'.  .  *=  T,  .  + 
i.J  i»j 


-At 


where  T'.  represents  the  temperature  at  an  advanced  tin*. 

J 

8.  Calculate  or  specify  1^  j  for  all  boundary  points  using  the 

applicable  boundary  conditions,  ..  •  -  previous  i.o  . 

9.  Calculate  using  Equation  (6-12)  and  the  T'^  j  just 

obtained  and  then  Z',  ,  for  all  interior  points. 

1  »  J 

10.  Calculate  S'j^j  using  Equation  (6-13)  for  all  Interior  points. 
A  Gauss-Seidel  point  iterative  procedure  is  used.  The  scanning  is 
performed  column  by  column  (i.e.,  from  j  =  2  to  jmax-1  for  i  «  2  to 
imax-1).  In  the  calculations  the  freshly  calculated  values  of  ,  are 
always  used.  The  over-relaxation  factor  W  is  calculated  from  the 
following  formula  (51,  Section  7.17): 


W 


1  +VTTX? 


,  vbsre  a-  f  +  ”  jS^l)  • 


The  iteration  is  terminated  as  soon  as  the  convergence  criterion: 


n+1  n 

Si,j  ’  Si.j 


’i.J 


&  EPS  for  every  (i,j) 


is  satisfied.  Here,  the  superscript  n  means  the  n-th  iteration  and  the 
value  of  EPS  used  for  all  calculations  was  0.001. 

11.  Calculate  Z'^j  and  S ' ^  j  for  all  boundary  points  using  the 
applicable  boundary  conditions  discussed  in  Chapter  VI. 

Hereafter  all  primed  quantities  refer  to  those  at  an  advanced 
time. 
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12.  Calculate  u'j  i  and  v'.  ,  for  all  interior  points  using 

>  J  i ,  j 

Equation  (6-14)  and  the  stream  function  just  calculated. 

13.  Determine  the  maximum  values  of  u'  and  v'  and  label  them  as 
uniax  and  vmax*  respectively. 

14.  Return  to  step  5  and  repeat  the  complete  cycle.  The  complete 
calculation  procedure  is  terminated  when  the  following  convergence 
criterion  is  met: 


Numerical  results  s 


Z’ 


Z . 


Zi,j 


£  10"4  for  all  (i,j). 


;how  that  when  this  criterion  is  met,  the  maximum 
error  calculated  in  step  10  is  usually  smaller  than  10“^,  and  all  the 
other  dependent  variables  (u,v,T,S,Z)  show  no  appreciable  changes. 

15.  The  output  consists  of:  partial  printout  for  intermediate 
results  during  the  calculation,  card  output  when  either  the  solution  is 
converged  or  the  execution  time  limit  is  exceeded,  and  complete  printout 
for  the  numerical  results  of  all  the  dependent  variables. 

The  computing  time  depends  on  the  combination  of  the  array  size  and 
the  source  size.  For  the  cases  studied,  the  computer  execution  time 
required  for  obtaining  a  converged  solution  increases  with  Gr  and  the 
source  size.  A  summary  of  the  execution  time  for  the  cases  run  are 


given  below: 
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No. 

Grid  System* 

Time,  min. 

Remarks 

Gr^  * 

JLO^ 

1 

13x15,1-4 

3 

2 

17x25,1-7 

11 

Continued  from  solution  of 
No.  1. 

3 

17x13,1-5 

6 

4 

17x17,1-5 

8 

5 

17x21,1-5 

10 

6 

17x25,1-5 

12 

7 

15x17,1*5 

5 

S 

29x33,1-9 

32 

Continued  from  solution  of 
No.  7. 

GrT  s 

105 

9 

15x17,1-5 

8 

Initialized  with  solution  of 
No.  7. 

10 

17x33,1-9 

14 

Continued  from  solution 

No.  9. 

11 

21x33,1-17 

17 

Continued  from  solution  of 

No.  10. 


•k 

See  p.  for  explanation. 


Appendix  G 


COMPUTER  PROGRAM  LISTING  -  FLOW  FIELD  CALCULATION 

This  computer  program  was  written  in  FORTRAN  IV  computer  language. 
Five  primary  variables  (the  dimensionless  indraft  U,  updraft  V,  tempera¬ 
ture  T,  stream  function  S,  and  vorticity  function  Z)  are  computed  for  an 
array  of  mesh  points  of  IMAX  x  JMAX.  There  are  four  major  parts  in  this 
computer  program:  reading  input  data,  calculating  all  the  repeated 
constants,  main  body  of  the  calculations,  and  output. 

The  input  listing  is  described  below: 

IMAX,  JMAX  maximum  values  of  I  and  J 

11,  J1  IMAX  -  1,  JMAX  -  1 

12,  J2  IMAX  -  2,  JMAX  -  2 

K  a  register  of  the  time  steps  advanced 

LL  total  number  of  iterations  performed  up  to  the  K-th  time 

step 

M  size  of  source  radius,  D/2  =  (M  -  1)  Ax 

MM  =  M  +  1 

N  frequency  of  intermediate  printout,  i.e.,  one  printout 

after  every  other  N  time  steps 
IM  number  of  columns  per  line  in  the  printout 

maximum  allowable  number  of  time  steps  to  be  advanced,  a 
limit  to  the  computer  time  that  may  be  used 


KM 
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LM 

II,  JJ 


GR 

PR 

EPS 

TIM 

TI 

XY 

UMAX,  UMAX 

CARDS 

READS 

W 

TEXC 

DIVIDE 


maximum  allowable  number  of  iterations  for  solving  the 
stream  function 

array  size  of  the  card  output  for  the  presently  converged 
solution  (This  array  will  be  used  as  the  input  for  finer 
mesh  spacing  continued  calculations.) 

Grashof  number 
Prandtl  number 

convergence  criterion  in  the  iteration  procedure  for 
solving  the  stream  function,  usually  set  at  0.001 

*  0.0,  for  initiating  the  dimensionless  time 

a  factor  for  adjusting  the  size  of  the  time  imJuunnil  to  Jbe 
advanced,  usually  set  at  1.0 

the  ratio  Ay/Ax,  1.0  was  used  for  the  present  calculations 
the  maximum  absolute  value  of  the  component  velocities  u 
and  v 

i  0.0  means  card  output  is  desired 

#  0.0  means  card  output  of  previous  calculations  is  to  be 
read  in  addition  to  the  regular  input  of  two  cards 
over-relaxation  parameter 

maximum  allowable  time  for  execution,  a  time  control  device 
^  0.0  mean  the  data  input  from  previous  calculations  is  to 
be  interpolated  to  obtain  the  values  at  the  new  additional 
points  for  a  subdivided  mesh  spacing  system 


skipping  instruction  for  boundary  value  specif ications : 
SKIP  >0.0  means  skipping  S(IMAX,  J),  T(IMAX,  J) ;  SKIP  * 

2.0  means  skipping  T(IMAXf  J),  S(IMAX,  J),  Z(IMAX,  J), 

T(I,  JMAX) 

convergence  criterion  (usually  set  at  0.0001),  if  met  the 

complete  calculation  is  terminated  and  the  complete  results 
are  printed  out 
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I5JMENSJCN  U ( 43 .43  )  ,  VJ43.43),  TI43.43),  Z<43,43>.  S(43,43).  DTOT ( 4 
13.431  .  DZDT(43.43).CA(43I .CBC43 > .CCI 43 ) ,CD( 43 > .CE < 43 ) ,CF < 43 ) »CG( 43 1 
l),CH(43).CH43).CJI43) 

201  FORMAT  </5X,  26HSTREAM  FUNCTIONS  —  SU.J)) 

202  FORMAT  {/5X.  2SHRA0IAL  VELOCITY  —  Utl.JM 

203  FORMAT  (//5X.  27HVERT ICAL  VELOCITY  —  V(I.JI) 

204  FORMAT  I//5X,  21HTEMPFRATURF  —  TII,J|) 

203  FORMAT  <//5X,  1 4HV0RT I CITY  —  7(T,J>) 

20fi  FORMAT  (1H1) 

207  FORMAT  (4X.4MK  -  .I4.4X.4HL  ■  . 1 3 . 5X , 7HT IME  »  , FI  2. 5 ,5X ,5HDT  •  , 
2E10.3.5X.5H0X  «  , F6. 3 ,5X  . 5HDY  *  .F6.3.5X.6HFPS  »  .F6.3.5X, 

35HTI  *  ,F6.3) 

208  FORMAT  (1H1.//.78H  NUMERICAL  SOLUTION  OF  NATURAL  CONVECTION  NEAR  A 
2  CIRCULAR  BOUNDARY  HEAT  SOURCE »4X , 4HGR  ».  F 10.3 »3X .4HPR  *,  F4.1, 

33X.7HGRID  «  .I2.1HX.I2I 
200  FORMAT  (33X.1PRF11.3) 

211  FORMAT  </) 

212  FORMAT  ( 16I5/E10.3.14F5.2) 

213  FORMAT  (6H  K  *  .I4.4X.4HL  «  . 1 3 .4X ,6HT IME  * .El 1 .4 ,4X »4HDT  *.Ell,4 
2 . 4X . 7H0MAX  «  .E11.4.4X. 7HC0NV  «  ,E1 1 .4.4X.5HLL  ■  ,18) 

216  FORMAT  (1P10E13.4) 

217  FORMAT  (13X.1P7F13.4) 

21<J  FORMAT  I1P8F10.3) 

270  FORMAT  (1PI2FH.X) 

221  FORMAT  (3X.1P7E11.41 

222  FORMAT  (2110,  1P4E15.7) 

C 

C  READ  DATA  INPUT  -  FOR  CALCULATIONS  AND  CONTROLS 

101  READ  (5,212'  IMAX.JMAX.il .J1.I2.J2.K.LL  ,M ,MM  ,N, I M ,KM ,  LM.II.JJ 
2. GR.PR.EPS, TIM  , T I ,X Y .UMAX , VMAX .CARDS .READS »W .TEXC .0 1 V  IDF ,SK I P  ,C0 
TCHK 1  ■  TIME(l.O) 

WRITE  (6.212)IMAX,JMAy,Il,Jl,I?,J?,X,LL,M,MM,N,IM,KM,  LM.II.JJ 
2. GR.PR.EPS, TIM  »T  I  ,XY .UMAX ,VMAX .CARDS .READS »W ,TFXC ,DI V I DF ,S< I P,CO 
CALCULATE  REPEATED  CONSTANT  TERMS  IN  MAIN  CALCULATIONS 
XM  .  M-l 
XN  *  N 
DX  »  0.5/XM 
DY  »  XY*DX 
DX2  *  DX*DX 
DY2  =  DY*DY 
C  =  1 . O/SORT ( GR ) 

Cl  «  C/PR 

C2  *  2.0*1 1.0/DX2+1.0/DY?) 

D  «  0.5*DX2/(l.0+DX2/DY2) 

DTI  *  C1*C2 
0XY2  =  7oO*DX2*DY2 
CI1-I1 
0  !  2*  T  2 

«  CT7*nx*Dv 
c 4  .  1  (  n*DX)**?/2.0 
1  DO  2  1. 2,  I  MAX 
CII  .  1-1 

COEFF  .  2. 0*C I  I *DX 
a i  «  0.5/rn 
*T  *  1.5/rri 
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ca<  n 
cbi  n 
ccf  n 

CD  (  I  ) 

cei  n 

CF(  I ) 

ccw  n 

chi  r  i 

CM  11 
CJ(  1 1 


( 1 .0-.AI )/DX? 
<1.0-AI)/0X2 
I  1.0+BD/DX2 
(  1.0-BI1/DX2 
c  r i i*Dxi#*2 
COEFF*OX 
COEFF*DY 
C I  I **2*DX Y2 
A  I / 1 DX*DY ) 
AI/DX2 


IF  (READS. NE. O.O) 
INITIAL  CONDITIONS 


GO  TO  107 


DO  6  I 
DO  6  J 
U( I  .J) 
V( I . J) 
T  (  I  ,  J ) 
2(1. J) 
S(  I  .J) 
DO  4  I 
T(  I  .1) 


*  1  •  I  MAX 

*  l.JMAX 

*  0.0 
*  0.0 
=  0.0 
.  0.0 
*  0.0 

■  1  ,M 
*  1.0 


GO  TO  82 

CARD  RFADING  FOR  CONTINUATION 

107  READ  IS. 222)  KK  »  LL  »  T I M .UMAX  »  VMAX  *GR 

IF  (D| VIDF.EO.O.O )  GO  TO  5 06 

C  INPUT  DATA  RFAOING  W5  "  SUBDIVISION  OF  MFSM  .SPACING  IS 
RFAD  (5.221)  1<U(I,J),V(I,J),TU,J),$U,J).ZH.J), 
2  J*  J ,  JMAX  .  2 1 


DFSIRFD 
1*1.1 MAX . 2 )  * 


10 


20 

21 


22 


00  20 
DO  20 
U( I .J) 

V  (  I  .  J  ) 

T  (  I  ,  J  ) 

S  (  I  *  J  ) 

2  (  I  .  J  ) 

DO  22 
DO  22 
U( I .J) 

V  (  I  » J  ) 

T  (  I  .  J  ) 

S  (  [  ,  J  ) 

2(1. J) 

T ( MM,  1  ) 

IF  (SKIP. NE. 2.0) 
PEAO  (5.219) 

READ  (5,219) 

RFAD  <t»219) 

RFAD  (5,219) 

PCAD  (5,210) 

RFAD  (5,210) 

GO  TO  ftp 


J«1 ,JMAX ,2 

I *2 » I  ma  x , 2 

<U( 1-1 ,J)*U( 1*1 .J)  1/2.0 
( v ( i-l .j ) ♦ v (  1*1 . J) i/2. n 
!T(  I-l  , J ) ♦ T (  1*1 ,J)l/2.0 
(5(1-1 ,J)*5(  I  +  1 ,Ji  )/2.0 
(  2  (  I-l  .J  )  +  Z (  1  +  1 . J)  ) /2.0 
J-2.JM5X  ,2 
1*1 .IMAX 

(  U ( I . J+l ) +U< I » J- 1 >1/2.0 
(  V(  I  .  J*1  )*V(  I  .  J-l  !  , /2.0 
(  T ( I ,J*1 )  +  T( I, J-l ))/2,0 
( S ( I , J+l )  +  5( 1 ,J-1  ) ) /2,o 
(2(1 ,J+1  )+2 ( I ,J-1  )  1/2.0 
0.0 

GO  TO  88 


( T ( JMAX.J) , 
( T ( I , JMA  X ) , 

( S< IM/X.J) , 

( S ( I , JMAX ) , 

(  Z( t MAX, J) , 
<Z(  I ,JMAX)  , 


J*2.16»2> 
1*2. IMAX, 2) 
J  =  2 ,4 , 2  ) 
1=2, IMAX, 2) 
J=  2 , 6 , 2  »  ) 

1  =  2, IMAX, 2) 
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88 


READ  (5,221) 
J=1  ,  'MAX  ) 

<  =  KK-i 

r.n  to  m2 


(("»'. J),V(I,J),T(I,J).S(I,J).Z(I.J),  1*1  ,  IMAX  1  , 


SIZE 


176 


r  MAIN  calcui at  tons 

C  TIME  TO  BE  ADVANCED  BASED  ON  STAB IL I  TV  CRITERION 
82  DT  ■  T I / ( UMAX /DA  +VMA  X/DY,OT ] ) 

TIM  *  TIM+OT 
BO  TCHK2  *  TlME(l.O) 

C  EXECUTION  TIME  CONTROL 

IF  (TCHK2-TCHK1.GT.TFXC  I  GO  TO  100 


K  * 

K  +  1 

IF 

(K.GT.KM)  GC  TO 

1  rr. 

ENERGY 

EQUATION 

9 

DO 

10  1*2.11 

DO 

10  J-2.J1 

IF 

(U(I.JI)  51.58,52 

51 

T1 

*  UC I »J)A( TC 1*1 . J> 

-  T ( I , J) )/DX 

GO 

TO  53 

58 

T1 

.  0,0 

GO 

TO  53 

52 

T 1 

«  U(i ,J)*(T( I , J)  - 

T ( 1-1 ,J) )/DX 

53 

IF 

(V(I.J))  54,59,55 

54 

T  2 

»  V(I.J)*(T<i,J*l) 

-  T( I ,J» )/DY 

CD 

TO  56 

59 

T2 

«  0.0 

GO 

TO  56 

55 

T? 

*  V(I»J)*iT(I,J|  - 

T ( I ,J-1 ) )/DY 

56 

IF 

( IcFO.l 1  GO  TO 

91 

T3  «  CA(1)*T(!*1*J)aCB( I ) *T ( 1-1  » J  )-C2*T  ( I  » J )  +  ( T ( I »J*1 > ♦T ( I  .J-l ) )/ 

1DY2 

RO  TO  K 

91  T3  »  4,0*<T(2,J)-T(l.J)>/DX2'MT<l,J-*-l)-2.0*m.J>'»T(l,J-in/DY2 

10  DTDT(l.J)  .  -T1-T2«<1*T3 

11  DO  12  1-2,11 
oo  12  j  «  ?,ji 

12  T(I,J>  >  T (  I  • J )  ♦  OT*DTOT ( I , J) 

27  DO  28  !«MM,I1 

28  T(I,1)  *  (18.0«T(  I.2)-9.0*TU»3)>2.0*T(  t«4)  1/11.0  NO  FLUX 

C  S  c  I P  »  2  »  0  ,  VALUES  ON  BOTH  IMAX  AND  JMAX  ARE  FIXED.  SKIP-1,0,  FIXFD  ON  IMAX  ON 

IF  (SKIP. CO. 2.0)  GO  TO  *3 

17  DO  18  I-  2,11 

18  T ( I  .UMAX)  -  T(  I  » J1  ) 

43  DO  44  J-2.J1 

T1J  -  (18.0*T(2.J)-9.0*T(3,J>*2.0*T<4,Jn/11.0 
IF  <.NOT.<T1J.OT.O,0,AND.T1J.LT.1.0>  )  GO  TO,  45 

IF  ( T 1 J, L T , T ( 2 , J) 1  GO  TO  45 

T  (  1  »  J  )  .  T  1  J 
GO  TO  44 

45  T  (  1  ,  J  )  -  T  (  2  ,  J ) 

44  CONTINUE 
C  VORTJCITY  FOUATION 

D7MAX  ■  0,0 

13  DO  14  1-2,11 

DO  14  J-2.J1 

IF  IUIIiJI)  61.68.62 

61  Z1  -  U(I ,J)*(Z(I+1 .J)  - 
GO  TO  63 


z ( i . j) ) /ox 
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7'  * 

r,A  rn  ^ 

21  *  uu.Ji*i?u,ji  -  zu-hjM/nx 
63  I F  (V(f,j)j  64,60,66 
66  17  --  Vrf  ,j)#(Z(!,j+1  ,  _  ZU,Jn/DY 
60  TO  66 
6«  72  *  O.n 
60  TO  66 

f c  *J’  -  ?< r.j-i » >/r>y 

66  IF  U.FO.l)  r,0  TO  6? 

l5?2"  CC<  n#z't  +  1,J>*cc><>)*2(l-i,j)-C2*2tl,j)  +  (Zil,j+ij42(jf  j_ln/ 

24  =  ITfl  +  I.J)  -  T I  I - 1 , J) ) »C J ( t ) 

60  TO  43 

03  OZOTd.j)  ,  -2 1  -z  2+("*2  3  +  24 

2FRR04  *  APSIOZDT ( I, J) 22(t , j) , 

16  IF  (2FRR0R.6T.D2VAX)  D2MAX*2FRROR 

15  00  16  I  =  ?,n 

00  16  j  J  2,J1 

16  2(1, Jl  =  2(1, J)  +  DT#DZDT(I,j| 

66  oo  47 

63  3-pt 

30  |_  -  n 

11  l  =  L  ♦  1 

IF  (L.6T.LM)  AO  TO  151 

RMAX  =0.0 

r  :T;^rw,r'""'  bv  si,rrf”1'"  o.n..wl„,„oNS 

00  24  J  =  ?,J J 

OSSOR  =  W*DS+( 1.0-W)»S( I ,J| 

IF  <6BMS((,ji,.LT,o.if.15)  ro  T 

0  =  AP4(D430R/S( I , J ) ) 

IF  (O.RT.RWAX)  RMAX  r  0 
26  S  (  I  ,  J  )  i  OSSOR 

00  =  ARS(RMAX-1.0)/0T 

IF  (00.6T.FPF)  0(1  TO  yj 

151  OMAX  =  00 
LL  =  Ll+L 

IF  (  SK  I P . fjT . 0. 0 )  00  Tn  71 

?R  00  30  J  =  2,J1 

30  $ ( I  MAX , J )  =  S ( 1 1 , J ) 

31  00  3?  j  =  ? , I  MAX 

IF  (SKI P, p0. ? .0 )  00  TO  32 

S(  ?■  JMAX  )  s.  FIT  ,J1  )  . ' 

HI  J*AX)  =  ? (  I , J 1  ) 

*  (  S  (  1 , 3  )  -  8.0*S(T,2n/  6H(  r  i 
0  50|  VINO  F0°  VFi  Oft  tipj 

3  3  00  34  I  r  ? ,  I  ] 

"O  74  J  =  ?,ji 

UU.Jt  ,  ri(I)  *(S(I  * J* 1 ) -S (  I,J-)  j  ) 


180 


34  V(t.J)  *  + 

105  WRITE  (6.213)  K .  L»  TIM  ,  DT  ,  OMAX.  D7MAX .  LL 
IF  ( < «LT « 30 )  00  TO  35 

IF  (DZMAX.LT, CO)  GO  TO  100 

35  DO  36  J  •  2.J1 

VI J  *  (18.0*V(2.J)-0«0*V(3.J)*2»0*V(4»J))/11.0  •  3-DT 

36  IF  (VIJ.OT.O.O)  V(l.J)  *  VIJ 

C  PREPARATION  FOR  CALCULATIONS  FOR  THE  NEXT  TIME  STEP 
UMAX  «  0.0 
3«  1)0  40  t»  2,  II 
DO  *0  J  «  2.  J1 
UIJ  *  ARStUtt.J)) 

40  IF  <U1 J.GT.UMAX)  UMAX  .  UIJ 

VMAX  «  0.0 

41  DO  42  I  *  1,  II 
DO  42  J  «  2.  J1 
VIJ  «  A P S ( V ( I ,J) ) 

42  IF  (VIJ.fiT.VMAX)  VMAX  *  VIJ 

r 

C  OUTPUT 

112  IF  (FLOAT(K)/XN.NE.FLCAT(X/N))  GO  TO  82 

102  WR1TF  (6.208)  GR . PR. IMAX , JMAX 

WRITF  (6.707)  X.L.TIM  .DT.DX.DY .FPS.TI 

WRITF  (6.20?) 

WRITF  (6,27")  (dl(T.J).  T  *  1  ♦  I M )  .  J.),JMAX) 

WRITF  (6.203) 

WRITF  (6.220)  ((V(I.J).  1  =  1. IM).  J-l.JMAX) 

WRITE  (6.204) 

WRITE  (6.220)  ((T(I.J).  I-l.IM),  J.l.JMAX) 

WRITF  (6.201) 

WRITF  (6,220)  ((SU.Jl.  1  =  1,  IM>,  J.l.JMAX) 

WRITF  (6,205) 

WRITF  (6,220)  ((Z(!,J),  I.l.JM),  J.l.JMAX) 

81  GO  TO  82 

100  IF  (IMAX.EO.IM)  GO  TO  104 

WRITF  (6.208)  GR, PR, IMAX, JMAX 

WRITE  (6.207)  X.L.TIM  ,DT ,DX ,DY .EPS .T I 

WRITE  (6,202) 

IF  ( IMAX. GT. 20)  GO  TO  110 

WRITE  (6,216)  IIUU.J).  1  =  1,10),  J.l.JMAX) 

WRITE  (6,211) 

WRITF  (6,217)  ((U(I.J),  I.ll.IMAX),  J.l.JMAX) 

WRITE  (6,203) 

WRITE  (6.216)  ((V(I.J),  1=1,10),  J.l.JMAX) 

WRITF  (6.211) 

WRITE  (6,217)  KVII.J),  !»11,IMAX),  J-l.JMAX) 

WRITE  (6,204) 

WRITF  (6,216)  ( ( T ( I » J )  ,  1.1,10),  J.l.JMAX) 

WRITF  ((,,711) 

WRITF  (6,717)  ( ( T ( I  ,  J )  ,  I.ll.IMAX),  J.l.JMAX) 

WRITE  (6,201) 

WRITF  (6.216)  (IS(I.J),  1.1,10),  J.l.JMAX) 

WRITE  (6,211) 

WRITE  (6,2)7)  ((S(!»J).  I.ll.IMAX),  J.l.JMAX) 

WRITF  (6.205) 
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WPITF  (6.216)  t (7 (  I  » J) *  T«1.10),  J»1 » JMAX ) 

WPITF 

WPITF  f  A  *  2 1 7  >  UZ(I.J).  I»1 1  * !MAX  )  ♦  J-l.JMAX) 

C, 0  TO  ]04 

1 1  0  WP  I Tr  (6.270)  ((UU.J).  T  ■  I  *12).  J»1  *  JMAX  ) 

WPITF  (6.211) 

WPITF  (6.70P)  l(Un.J).  I  »1  3  *  T  VAX  )  *  J-  l.JVAX) 

WPITF  f  A  *  2^3  5 

WPITF  (6,720)  ( ( V  (  I  ♦  J  )  »  T  *  T  *  T  2  >  * 

WPITF  (6.211) 

WPITF  ( 6 ,  2rt<?  )  ((V(I.J).  T  «  1 3  ♦  T  VAX  1  •  J«1*JMAX) 

WPITF  (6.704) 

WPITF  (6,720)  I'Tll.J),  I«1.!21.  J«1,JMAX) 

W9ITF  (6,,’lU 

WPITF  (6.70°)  ((  ,(J.J),  T«13.!VAX).  J«1»JVAX) 

WPITF  (6.?n)) 

W®  T  T  F  (A,??")  (fG(T,J),  J«1,JMAX1 

WPITF  U,?H) 

WPITF  f  6  *  2^0  I  (  (  5  (  I  «  J )  «  I  *  1  3  ♦  I  VAX  )  »  J-l.JVAX) 

WPITF  16.205) 

WPITF  (6.220)  ((Zd.Jl.  1*1.12).  J« l.JVAX) 

WRITE  (6.211) 

WPITF  (6.700,  ((7.(1. J).  T  »1  3  *  I VAX  )  .  J-1.JMAX1 

ln4  IF  (CAonF.FO.Fi.n)  GO  TO  103 

WPITF  (7.222)  K.LL.TIV.UVAX.VVAX.G® 

IF  (07MAX.LT.rO)  GO  TO  !"« 

WPITF  (7.221)  ((U(I.J)»V(T,J).T(!.J).5(T.J).7(I.J)»  T  *1  .  T VAX )  . 

2  J*1.J«AX) 

GO  TO  103 

1 0B  WPITF  (7.271)  ((U(I.J).V(I,J).T(I.J).S(I.J).7(I.J)»I»],II).J«1,JJ) 
in'.  r.O  TO  )o, 
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13  ABSTRACT 


The  turbulent  free  convection  of  air  above  a  2-foot  diameter,  heated  horizontal 
plate  has  been  studied  experimentally  and  numerically.  The  mean  temperature  fields 
and  the  indraft  profiles  for  two  mean  plate  temperatures  were  measured  using  a  thermo¬ 
couple  and  a  constant  temperature  hot-wire  anemometer.  Also,  the  turbulence  and  mean 
velocity  were  measured  for  the  higher  plate  temperatun  using  the  hot-wire  method.  The 
flow  field  was  visualized  by  shadow  photograph  technique.  From  visualization  and 
measurements,  it  was  found  that  the  region  of  significant  deviation  from  ambient 
temperature  and  velocity  was  restricted  to  a  region  near  the  plate  centerline  (the 
primary  flow  region).  The  indraft  velocity  was  found  to  be  relatively  large  near  the 
ground  level  (within  approximately  1"  of  the  ground). 

The  major  temperature  drop  took  place  in  the  region  very  near  the  plate.  Within 
0.02"  of  the  plate  the  temperature  distribution  in  the  air  could  be  calculated  based 
on  conduction  only.  This  region  was  therefore,  called  the  "conduction  layer."  At  a 
given  mean  plate  temperature,  the  temperature  gradient  was  found  to  increase  with  the 
radius.  Data  obtained  from  heat-transfer  measurements  were  consistent  with  the  one- 
third  power  correlation  reported  in  the  literature. 

The  turbulence  in  the  flow  field  was  found  to  consist  of  low  frequency  and  high 
amplitude  fluctuations  (on  the  order  of  10  Hz  and  1  ft/sec).  Because  of  the  limita¬ 
tion  of  the  hot-wire  technique  for  large  turbulence  measurements,  flow  velocities 
could  not  be  deduced  directly  from  hot-wire  data.  To  remove  this  difficulty,  a 
numerical  data  simulation  scheme  has  been  developed  in  which  the  parameters  describing 
the  turbulence  flow  (r.m.s.  fluctuations  and  correlation  coefficients)  were  used  as 
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DD,FrJ473  ,PAGE  " 

S/N  010 1.807. 6801 


UNCLASSIFIED 

Set  unis  t  Id •  si{h  c*ti«*n 


Input.  By  interring  from  the  simulated  data  of  known  parameters,  experimental 
hot-wire  data  reduction  was  then  possible.  Data  reduction  model  was  validated 
by  numerical  experiments. 

The  eddy  diffusivity  in  the  region  away  from  the  conduction  layer  was 
estimated  based  on  temperature,  velocity  and  turbulence  data  using  two 
independent  methods.  The  agreement  was  good.  The  spatial  variations  of 
the  eddy  diffusivity  in  most  of  the  primary  flow  region  was  found  to  be 
gradual  with  rapid  drops  occurring  in  the  region  between  the  primary  flow 
and  the  cold  ambient. 

A  numerical  flow  calculation  was  made.  The  mathematical  formulation 
was  baaed  on  Boussinesq  approximations  using  a  constant  eddy  diffusivity 
model.  A  turbulent  Grashof  number  Gr>j>  (the  governing  parameter)  was  defined 
through  the  definition  of  a  characteristic  plate  temperature  rise  A#T,  the 
plate  mean  heat  flux  and  the  eddy  diffusivity.  Gr^  and  A@j  were  obtained 
based  on  the  best  fit  of  experimental  and  numerical  centerline  temperatures. 

By  the  specification  of  AAp  at  the  plate  surface,  the  effect  of  the 
intense  variation  of  eddy  diffusivity  in  the  conduction  layer  region  could 
be  avoided  in  the  numerical  calculations.  Numerical  results  based  on  a 
constant  eddy  diffusivity  model  were  obtained  and  compared  with  the  experi¬ 
mental  data.  Due  apparently  to  the  non-constancy  of  the  eddy  diffusivity, 
the  calculated  temperature  and  velocity  profiles  exhibit  less  constriction 
than  the  experimental  data.  Therefore  a  more  general  turbulent  transport 
model  will  be  required  to  provide  a  good  theoretical  description  of  the 
phenomena. 
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